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I. INTRODUCTION 



Stellar mass compact objects inspiraling into super- 
massive black holes (SMBHs) will be important astro- 
physical sources of gravitational waves (GWs) for future 
space-based detectors. Accurate GW templates for such 
extreme mass-ratio inspirals (EMRIs) require detailed 
knowledge of the motion of the source, so there has been 
a community effort to calculate EMRI trajectories. If 
we neglect the gravitational self-force of the small ob- 
ject, its orbit is a Kerr geodesic that, up to parameters 
specifying the initial position, is characterized by three 
constant orbital parameters: an energy E, an azimuthal 
angular momentum L z , and the Carter constant Q. De- 
termining the inspiral is tantamount to calculating how 
the self-force causes both the positional parameters and 
the orbital parameters to evolve in time. 

Despite ongoing efforts, direct evaluation of the self- 
force in the Kerr case is still not possible. Accordingly, 
there have been parallel efforts to approximate its ef- 
fects. The focus of this paper is the adiabatic approx- 
imation, which captures the slow secular evolution of 
E, L z , Q by solving a system of ordinary differential equa- 
tions (ODEs) of the form 

dE 

— =F E (E,L Z ,Q) (la) 

d Jg=T Lz {E,L z ,Q) (lb) 

d ^-=T Q {E,L z ,Q) . (Ic) 

For now, it suffices to know that the righthand sides 
(RHSs) of equations ([I]) are so costly to evaluate that 
these equations will have to be integrated using a nu- 
merical grid. More specifically, the EL Z Q velocity field 
will be pre-computed only on a dense mesh of points in 
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FIG. 1: Above is a heuristic depiction of two possible nu- 
merical grids that could be used to generate adiabatic ap- 
proximations (dashed curves) to true inspirals (solid curves) 
in the orbital parameter space. The dots represent a set of 
resonant grid points and the plus signs a set of non-resonant 
grid points. The resulting adiabatic curves are the same in 
either case but significantly less costly to produce with the 
resonant grid. A true inspiral may evolve in a way that is not 
well-approximated adiabatically as it approaches a low-order 
resonance, as on the left. That divergence, if it occurs, hap- 
pens regardless of whether the resonant point is used as part 
of the numerical grid. 



EL z Q-sp&ce. Real-time integration of ([I]) will then rely 
on derivative values interpolated off of that grid. 

In this paper, we advocate building such grids us- 
ing only points corresponding to resonant geodesies, for 
which the frequencies of the radial and polar motions 
are rationally related. As we will see, intermediate cal- 
culations that comprise the bulk of the computational 
expense can be recycled among several Fourier modes on 
resonant grid points but must be recomputed for every 
mode in the non-resonant case. We estimate that, com- 
pared to using non-resonant grid points, our prescription 
could reduce the computational cost of an EMRI grid by 
an order of magnitude or more. The resonant-grid pre- 
scription will also facilitate faster computation of GW 
snapshots from geodesic sources. 
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We represent our proposal schematically in Figure [T] 
First, the RHSs of equations are evaluated directly 
on a grid of either resonant points (dots) or non-resonant 
points (plus signs). At any other point, the RHS values 
can be interpolated from the values at the grid points. 
The adiabatic eqautions are continuous and smooth, 
so regardless of which grid is used, integrating them pro- 
duces the same adiabatic solutions (dashed curves). The 
only difference is that the resulting adiabatic curves cost 
significantly less to generate with the resonant grid. 

Ref. p] and more recently, Ref. [2] have noted that 
such adiabatic approximations may fail to capture impor- 
tant features of the true inspiral (solid curves) near low- 
order 1 resonances. Heuristically speaking, those authors 
argue that while an adiabatic solution may remain fairly 
faithful to an inspiral that steers clear of resonant points 
(lower right of the figure) , those approximations may fare 
much worse for an inspiral that transits near a resonant 
point (middle left of the figure). To pre-empt possible 
confusion, we remark that there is no inconsistency be- 
tween this observation and our proposal. The decision to 
include any particular EL Z Q point, resonant or not, in 
the numerical grid is unrelated to whether the resulting 
adiabatic curves will faithfully reflect EMRI motion near 
that same point. The ironic coincidence is that the points 
where the adiabatic approximation is most likely to fare 
poorly 2 are also the optimal grid points for generating 
adiabatic curves. 



II. RESONANT KERR ORBITS 

The paramount role of resonant orbits was the central 
theme of an earlier series of papers [8-10J. (We use the 
terms "resonant", "closed", and "periodic" interchange- 
ably.) A spectrum of closed orbits, which manifests as a 
spectrum of multi-leaf clovers, entirely structures black 
hole dynamics. Although completely closed orbits must 
return to their initial values 3 of (r, 9, ip) simultaneously, 
only the r-9 periodicity detailed in Ref. [S] is relevant to 
the present work. The rational number associated with 
the r-9 frequencies determines the multi-leaf clover ge- 
ometry. What's more, that rational obediently stacks 
in energy monotonically: lower rationals correspond to 
lower energies than do higher rationals. Such orbits con- 
stitute a measure zero set but are nonetheless dense in 
the phase space. Every non-resonant orbit is arbitrarily 
close to some resonant orbit. 

We first consider resonant orbits in physical space and 
then again in phase space. Since we will be concerned 
with functions that do not depend explicitly on either 
the azimuthal angle or on coordinate time, it will suffice 
for us to restrict attention to geodesic motion in two co- 
ordinates (r, 9) in physical space and to the projection 
of the motion into a AD submanifold of the phase space 
spanned by (r, 9) and their conjugate momenta. 



A. Resonant Orbits in Physical Space 



The rest of this paper is organized as follows. In Sec- 
tion [TIJ we review some relevant features of resonant Kerr 
orbits in both physical space and phase space. In Sec- 
tion [nij we summarize how one arrives at the adiabatic 
equations of motion and clarify why an averaging per- 
scription required to derive those equations must be a 
torus average rather than a time average, an issue that 
has raised some debate in the literature [IJ[3HZ]- Partly to 
help make the averaging argument and partly because we 
will focus on frequency-domain approaches to solving the 



adiabatic equations, Section III also provides some nec- 
essary mathematical background on Fourier analysis in 
both the non-resonant and resonant cases. With a clear 



view of the adiabatic program now in hand, Section IV 



presents the main result of the paper, namely a concrete 
prescription for computational savings that frequency- 
domain EMRI codes can leverage by using resonances. 
Finally, Section [V] speculates about how a more unortho- 
dox use of resonances could offer additional efficiencies 
provided it can be practically implemented. 



The black hole is completely characterized by its 
mass M and spin a. The geodesic of the lighter com- 
panion is characterized by four dimensionless constants 
H, E, L Z ,Q. In Boyer-Lindquist coordinates and dimen- 
sionless units, which is equivalent to setting M = [i = 1, 
the radial and polar Kerr equations of motion can be 
written as 



= ±y/R(r) 

= ±V®¥) 



(2) 
(3) 



where 



R(r) = -(l-E 2 



0(9) = 
and 



2r 3 (4) 
[a 2 (1 - E 2 ) + L 2 ] r 2 + 2 (aE - L z f r - QA 

^ i - i ( 5 ) 



Q- COS 2 (9)!a 2 (l-E 2 ) + ^f-} 
I sin 9 J 



r 2 - 2r + a 2 . (6) 



A resonance is low-order if the numerator and denominator of 
the rational frequency ratio are both small integers. 
2 To balance the argument, Ref. [I] also offers plausible reasons 
why the adiabatic approximation may still be valid near reso- 
nances. 3 ip need only return to its initial position mod 2ir. 
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An overdot denotes differentiation with respect to 
Mino time 11 , A, related to proper time r by 



dr 
dX 



= £ = r 2 + a 2 cos 2 ( 



(7) 



The advantage of using Mino time is that the r and 8 
equations of motion decouple and that each is a function 
only of one coordinate. To make connections with obser- 
vations, we will often care about how certain quantities 
evolve with respect to coordinate time t. However, coor- 
dinate time turns out to be mathematically cumbersome, 
so throughout this paper, we perform all intermediate 
calculations related to such quantities by first changing 
variables to Mino time. 

Solving equations Q and ^ for the radial and polar 
turning points, we find that the radial coordinate varies 
between a periastron r p and an apastron r a and that the 
polar coordinate similarly varies between some minimum 
value Omin and maximum value 9 max — n — 9 min . All 
turning points depend only on the constants E,L Z ,Q. 
We introduce the simplifying notation 



£ = (E,L Z ,Q) 



(8) 



for those 3 orbital parameters and reserve the symbol £ 
to refer to any one of E,L Z ,Q individually. 

The radial and polar coordinates are each periodic with 
respective Mino periods 



A r = 2 
A fl = 4 



dr 



^ 2 d9 



le min VW) 
and corresponding frequencies 

2tt 



A,. 
2k 

a7 



(9a) 
(9b) 

(10a) 
(10b) 



The radial and polar velocities are also periodic with the 
same corresponding periods and frequencies. If the fre- 
quency ratio 



l + q r e 



(11) 



is a rational number, an r-9 projection of the resulting 
orbit closes after a finite time. 

Note from equations Q , ^ and ([9| that the frequen- 
cies and q r g depend only on the constants £. q r g is also 
a topological invariant and thus coordinate independent. 
A q r g for which the relatively prime numerator and de- 
nominator are both low-valued integers will be referred 
to as "low-order" . We arbitrarily call low-order resonant 
orbits those for which the numerator and denominator of 
q r g are each less than 10. 
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FIG. 2: Top: A low-order periodic orbit with q r e = i, a = 0.9, 
E = 0.954788, L z = 2.65115 and Q = 0.944969, projected 
on the r-9 plane (we plot r-cosO to make the figure more 
viewable) with initial conditions ro = r a = 17.8148 and #o = 
1.22079. Bottom: A non-resonant orbit with q r g ~ 
^jTjiur v a = 0-9, E = 0.956, L z = 2.65115 and Q = 0.944969, 
with initial conditions ro = r a — 18.4568 and 9o = Omin = 
1.22076. 



Projections of periodic orbits into the r-9 plane pro- 
duce Lissajous figures. The top panel of Figure [2] shows 
the Lissajous figure of a periodic orbit with a low-order 
q r e, while the bottom panel shows the analogous projec- 
tion of a neighboring orbit with an irrational q r g. 

The figures produced by projecting into the r-9 plane 
are less topologically insightful than the figures in an or- 
bital plane, loosely defined in the Kerr system as the 
plane perpendicular to the orbital angular momentum 
[S]. In the orbital plane, the rational q r g has power- 
ful topological information and can be interpreted as 
q r g = w + v/z, where the integer w represents the num- 
ber of nearly circular whirls near periastron, the integer 
z is the number of elliptical leaves in the multi-leaf clover 
pattern, and the integer v is the order in which the leaves 
are hit [El [9]. To illustrate, the same two orbits of Figure 
[2] are plotted in the orbital plane in Figure [3j The orbit 
in the top panel of Figure [3] has q r g — 1/2 and therefore 
corresponds to a 2-leaf clover, as is now evident. The bot- 
tom panel non-resonant orbit is close to a resonant orbit 
with q r g « 25o'ooo wmcn would correspond to a 250, 000- 
leaf clover that skips 125,857 leaves in the pattern each 

so that the orbit 



time. Notice that 250000 = \ 



sr,7 



250,000 

is really a tight precession of the 2-leaf clover through an 
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FIG. 3: Top: The same orbits from Figure [2j projected into 
the orbital plane. 



angle of 250 ooo J ^ n ~ 0-02154 radians per radial cycle. 

While they do fix the turning points, the constants £ 
do not uniquely fix the orbit [12]. An orbit that hits 
apastron at 9 m i n is not identical to the orbit that hits 
apastron at a different value of 9, as shown in Fig. [4] 
Since q r g depends only on constants, a q r g = 1/2 orbit 
is always a 2-leaf clover in the orbital plane [9]. How- 
ever, orbits with different r-9 initial conditions (r ,#o) 
are rotated relative to each other in the orbital plane. 

As Fig. [4] shows, the resulting orbits are genuinely dis- 
tinct in 3D. Presumably, they could have distinct grav- 
itational wave emissions. Interestingly, though, perihe- 
lion precession happens on a faster time scale than plane 
precession. It is therefore reasonable to suspect that all 
orbits with the same q r g generate similar waveforms and 
that the different plane precessions induce modest dif- 



ferences in the modulations of the amplitude [T3]- We 
remain agnostic on the relative importance of r-9 initial 
values on the waveform generated and instead focus on 
efficient calculation of adiabatic inspirals. 

B. Resonant Orbits in Phase Space — Phase Space 

Tori 

In the 4D space spanned by (r, p r , 9,pg), all orbits (res- 
onant or non-resonant) lie on 2D tori that can be con- 
structed as the Cartesian product of two closed curves. 
We obtain one of those closed curves if we project an orbit 
into the r-p r plane. The area of the curve is the familiar 
action J r used in action-angle coordinates. Analogously, 
the projection of the same orbit into the 9-pg plane yields 
another closed curve with area Jg. We now consider that 
pair of curves as a locus of points on a 2D surface with 
the topology of the 2-torus S 1 xS 1 = T 2 . Every set of or- 
bital parameters £ defines one such torus that we denote 
T|. 

The use of Mino time as an evolution parameter fur- 
nishes one (but certainly not the only) coordinate sys- 
tem for T|», according to the following construction. As 
already mentioned, the motions r(A) and 9(X) are each 
individually periodic in Mino time, with periods A r and 
Ag (and frequencies fl r and fig), respectively. Scaling 
the evolution parameter A on each of the r-p r and 9- 
pg curves by fl r and fig, respectively, leads to a natural 
definition of angle variables Xr = £l r X and xe = fleX. 
We choose a specific trajectory (r(X),p r (X),9(X),pg(X)) 
in order to assign \r and xe values, respectively, along 
the r-p r and 9-pg curves, but the trajectory is only a de- 
vice that we discard once the torus coordinate system is 
in place. The points at and 2-7T in each of \r , X6 are 
identified, so the torus can be represented as a 27r-by-27r 
square with opposite sides identified as in Fig. [5} We will 
make a simplifying choice that (r a ,# m i n ) corresponds to 
the origin 4 of the torus. Then, a reflection in the line 
Xr = 7r corresponds to keeping r fixed and reversing the 
sign of p r , and analagously for reflections in xe — ti\ 
Note that each quadrant of the toroidal square therefore 
contains the same (r, 9) pairs but with all possible sign 
combinations for the momenta (++, H — , — h, ). 

Note that each Xr corresponds to an ordered pair 5 
(r, p r ) and each xe corresponds to an ordered pair (9,pg). 



4 Many references, including [Tl 131 [TT1 114| . instead tacitly choose 
the point (r = r p ,p r = 0,9 = m in,Pe = 0) as the origin of 
the torus coordinates. We say "tacitly" because they refer to 
the individual orbit with those initial conditions as a fiducial 
geodesic to use in their analyses. Another interpretation of that 
choice is that they are working not with one geodesic but with 
one torus and that they have instead chosen a fiducial origin for 
a Xr~Xe coordinate system on that torus. 

5 Some references describe the mapping of functions of the form 
F(r, 0) to corresponding functions F(x r , Xe )• In fact, no function 
that enters an adiabatic EMRI calculation depends on r and 6 
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FIG. 4: The above figures are all q r g = | orbits with a = 0.9 and orbital parameters L z = 2.65115, Q = 0.944969 and 
E — 0.954788. The three figures in each row have the same initial coordinates. The column on the left shows an r-cos 9 
projection of the orbit, the middle column is a projection in the orbital plane, and the right column is the 3D orbit. All three 
rows have ro = r a = 17.8148. The first row has 6q = 9 min = 1.22079 and is the same orbit shown in the top panel of Figs. [2] 
and|3j the middle row has #o = 1.39579, and the bottom row has 9o = \. 



We discuss alternative coordinate systems for T|. in Ap- 
pendix [X] and elsewhere in this article but will use only 
the (XriXs) coordinates for calculations. 

On the compact (\r,Xe) square defined above, 



alone. The notation F(r, 9) in those references is used because, 
once restricted to a torus, the value of each coordinate determines 
its conjugate momentum up to a sign. Still, the values of those 
signs affect the value of the function. We believe a notation such 
as F(r,p r , 6,pg) for these pre-torus phase space functions is more 
appropriate. 



geodesic trajectories are lines of slope Q,g/Vl r = 1 + q r g. 
With respect to Mino time, those orbits are given para- 
metrically by 



X r (A) = Q. r X + Xr„ 



(12a) 
(12b) 



Two different initial positions xo and x'o produce dis- 
tinct orbits unless there exist real numbers x and y that 
simultaneously satisfy 



y-xe a 

X - Xr„ 



(13) 



(> 




FIG. 5: The above picture shows a resonant torus mapped to 
a square with the path of two resonant orbits traced out. The 
solid line shows the path of a resonant orbit with Xr = Xe = 
and the orbit traced out by the dotted line has Xr — an d 
Xe — 0.7894. The resonant torus and both resonant orbits 
have ^ = l + q r9 = B = §. 



and 



Xr = X m0d 2lT 

Xe = V m °d 27r 



(14) 



If these conditions are met, then the two different initial 
positions produce time-translated versions of the same 
orbit. 

When ilg/fl r is irrational, we will call both the torus 
and any orbits on that torus non-resonant. Orbits on 
non-resonant tori never close and instead sample the en- 
tire torus ergodically: an orbit starting from any ini- 
tial condition will pass arbitrarily close to every point 
in the torus after some finite (but possibly very long) 
time. Therefore, non-resonant orbits with different xo = 
(Xr ) XSo ) are arbitrarily close to time translations of ev- 
ery other non- resonant orbit with the same xo- We will 
alternately refer to such orbits as aperiodic or bipcriodic. 

When the frequency ratio Q,g/VL r is a rational number 
£ , we will call both the underlying torus and orbits on 
that torus resonant. Orbits that live on resonant tori 
inherit the rational frequency ratio and thus always trace 
out closed curves. Since no single resonant orbit ergodi- 
cally fills the torus, even after infinite time, two resonant 
orbits with the same £ but different {xr , Xo ) are n °t 
necessarily time translations of each other. The set of 
all resonant orbits with the same £ does fill the entire 



torus. Because they return to their initial conditions af- 
ter a finite time, we will alternately refer to these orbits 
as periodic or singly periodic. 

Figure [5] shows two resonant orbits on the resonant 
torus defined by E = 0.954788, L z = 2.65115, Q = 
0.944969. Each is thus a q r g = \ orbit, or 2-leaf clover 
in the orbital plane. These are the same two orbits il- 
lustrated in physical space in the top two rows of Fig- 
ure [4j The two orbits are distinguished by their initial 
position xo on the torus. The solid line orbit, which 
starts at Xr a — 0, Xe Q = 0, corresponds to the physical 
space orbit with initial conditions = r a = 17.81477 
and 6» = 6> min = 1.220793. The dotted line orbit with 
intitial conditions Xr = and xo Q — 0.7854 has phys- 
ical space intial conditions of r$ — r a = 17.81477 and 
do = 1.39579. Notice that any two adjacent line seg- 
ments belonging to a single orbit are separated in Xr by 
^ and in xe by =^ but are not traced out sequentially 
for general -. 

In the same way that the rational numbers have zero 
measure on the line, the set of resonant tori has zero 
measure in the 4D phase space. To date, most of the 
literature on the adiabatic EMRI problem has ignored 
resonant geodesies precisely for this reason. Nevertheless, 
as we will see, the judicious exploitation of this measure 
zero set leads to significant computational efficiencies in 
adiabatic EMRI calculations. 



III. AVERAGING IN THE ADIABATIC 
APPROXIMATION 

Given the background on geodesic dynamics, we now 
turn to the adiabatic approximation of EMRIs, an ap- 
proximation that has seen substantial debate in the lit- 
erature. As we elucidate below, most of that debate con- 
flates the question of what kind of averaging procedure to 
use in the equations of motion ([I]) with other related but 
logically independent questions about the adiabatic ap- 
proximation. In this section, we clarify why phase space 
averaging (as opposed to time averaging) is the correct 
averaging procedure. We also establish the results we will 
need in Section [IV| to exploit resonant orbits for compu- 
tational savings. 



A. The adiabatic equations of motion 

Let X denote the Boyer-Lindquist coordinates of the 
inspiraling object along with its canonical radial and po- 
lar momenta. In the absence of radiation reaction, the 
equations of motion are 



dX 
Itt 
d£ 



G(X,£) 



(15a) 
(15b) 
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where the RHSs G of the positional equations are some 
form of the equations for geodesic motion, e.g. Hamil- 
ton's equations for free-particle Kerr motion. Radiation 
reaction adds to the RHSs new functions 



dX 
~dt 
d£ 

~dt 



= G(X,£) + F(X,£) 
= + f(X,S) 



(16a) 
(16b) 



that are determined by the full gravitational self-force on 
the particle. Those unknown functions can be expanded 
in a perturbation series in powers of a natural small pa- 
rameter: the system's mass ratio ee/i/M«1. Further- 
more, at each order in e, the functions above decompose 
into a sum of dissipative and conservative pieces: 



F 



F 



(i) 



diss 



F {2) + 

diss 1 cons 



.#(2) 



m 

I cons 



■OP) 

m + H2) 

J diss 1 J com 



O £ 3 



(17a) 



(17b) 



See [7] and references therein for a fuller account. 

We expect a natural separation of timescales in this 
system. The "fast" positional variables X will change 
substantially on a short timescale equal to an orbital pe- 
riod T or b ~ M, while the "slow" orbital parameters £ 
only change substantially on the much longer timescale 
Trad ~ Mje. Due to the coupling of the X and £ 
equations, both X(t) and £(t) should exhibit oscillations 
around a secularly trending central value, but the oscilla- 
tions in £ should be 0(e). In such a system, a first-order 
averaging procedure seeks an approximate and hopefully 
more tractable set of equations for the slow variables from 
which the dependence on the fast variables (and thus the 
source of the small oscillations) is removed [15l - fT7] . 

Averaging must theref ore decouple the £ equations 
from the X equations in ( 16 ) in order to isolate the sec- 



ular trend in the former. 6 One can even adopt the point 
of view that the desideratum of a preliminary averag- 
ing procedure is to decouple the equations for the slow 
and the fast variables from each other as much as pos- 
sible. We represent those averaged, decoupled equations 
for the orbital parameters (equivalent to equation (fl])) as 



d£ 



secular 



df 




= (f(X,£ 



= T[£ 



(18) 



6 Note that the converse is not possible, since the fast variables are 
coupled to the slow ones at zeroth order, where the slow variables 
appear as constant parameters. 



Throughout this paper, we will represent averages of all 
sorts with angle brackets () and use subscripts on the 
brackets to denote the type of average implied. Note 
that in (18), we have used (• • •) to denote an average 



without yet specifiying which variables that average is to 
be taken over. 



B. Flux balance and its relationship to averaging 

Although we now have the general form of the adi- 
abatic equations, we cannot write them explicitly be- 
cause we still do not know how to evaluate the self-force. 
Mino showed [HE] that, under the assumption of non- 
resonance, the infinite time-averaged values 7 of the func- 
tions Te and Fl z equal the sum of the infinite time- 
averaged fluxes of the corresponding orbital parameters 
at radial infinity and the central black hole horizon in 
GWs emitted by the system 8 . While there is no con- 
served Q-current to associate with a GW Q-flux, Mino 
likewise showed that there are analogous infinite time- 
averaged quantities at infinity and the horizon that sum 
to the infinite time-averaged value of Tq. Though not 
strictly physically accurate, we will henceforth refer to 
those quantities as fluxes of Q for ease. Subsequent work 
[31[TH1[IS] has led to explicit formulae for these Q-fluxes. 

Fortunately, we do know how to calculate the afore- 
mentioned time-averaged fluxes at infinity and the hori- 
zon via the computationally mature Teukolsky formal- 
ism. Various Teukolsky-based (TB) codes can compute 
the required fluxes from equatorial orbits [20] ■ spherical 



orbits J (constant 



], and now generic orbits of arbi- 



trary inclination and eccentricity [TH [T^l HI] . 

These developments have led to the following two-stage 
implementation of the adiabatic approximation. The first 
stage, usually called the radiative approximation 10 , keeps 
only the lowest-order contributions from the dissipative 
self-force (since the conservative contributions will aver- 
age to zero). The second stage, called the flux-balance 
method, uses the time-averaged nonlocal fluxes (com- 



7 The quantities of physical interest are averages over coordinate 
time t. 

8 Mino's proof suggests that this equivalence is only true for 
non-resonant geodesies and possibly a small subset of resonant 
geodesies. 

9 Many authors refer to orbits of constant radial coordinate r as 
"circular" even when they are nonequatorial. We prefer the term 
"spherical" for such orbits (as used in 1211 ) and reserve the term 
circular for constant r equatorial orbits. 

10 There is some dispute about whether the neglected secular effects 
of the averaged conservative piece of the self- force manifest them- 
selves at the same order in the small expansion parameter fi/M 
as the dissipative pieces [3HZ]- That dispute does not concern 
us here. Whatever its limitations, the radiative approximation, 
is here to stay for at least the foreseeable future, if for no other 
reason than that it is the only relativistically correct approxima- 
tion to the inspiral motion accessible to numerical calculation in 
the status quo. 
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putable) as proxies for the averaged local contributions 
of the dissipative self-force (not currently computable). 
The RHSs of equation ( 18 ) end up with nonlocal fluxes 



inside the brackets and an interpretation of those brack- 
ets as infinite time-averages. 

There is a problem, however, with this perscription, 
which intertwines two logically distinct facets of the adi- 
abatic approximation to EMRIs: 

1. Is it mathematically appropriate to interpret the 



angle brackets in equation (18) as a time-average, 



or is some other sort of average required? 
2. Given the answer to 1, can we evaluate the RHS 



of (18), either directly or by finding a numerically 
equivalent proxy? 

After all, the fact that we can compute a time-averaged 
proxy does not imply that we should be time averaging 
in the first place. 



The form of equation ( 18 ) suggests two ways to average 



the RHS in order to remove the dependence on the po- 
sitional variables: for fixed £, we can either phase space 
average over the torus, or we can evaluate the RHS along 
a specific orbit on that torus and then average over time. 
In Section |III E[ we offer a definitive argument in favor 
of torus averaging instead of time averaging. 

To arrive at that conclusion, we must first distinguish 
between torus functions and time functions. Torus func- 
tions assign a value to every point on a phase space torus, 
while time functions assign a value to points along an 
individual orbit that are labeled by the value of an evo- 
lution parameter (i.e. a time variable). Our conclusions 
about adiabatic averaging will be based on differences 
in how Fourier analysis is done on these two domains 
— a 2-dimensional compact position space for the torus- 
functions and a 1-dimensional noncompact time axis for 
the time-functions. Moreover, numerically accurate flux 
calculations require frequency-domain TB codes that sep- 
arately compute fluxes from individual Fourier modes, 
and the aforementioned different domains also impact the 
details of the modewise flux calculation. 

Before delving into those details, we must mention an 
important point. Average values, whether in the torus 
or time sense, are coordinate dependent, and in certain 
applications it matters which coordinates the average is 
taken over. The angle brackets in equation (18), for in- 



stance, will turn out to denote a torus average not over x 
but over a different set of torus coordinates 7 = (7,-, 70) 
described in Appendix [XJ However, torus averages with 
respect to x are much easier to compute than those over 
7, in much the same way as Mino time averages are easier 
to compute than are coordinate time averages. Luckily, 
for every torus function U (7) and every time function 
u(t), we can always construct different functions V(x) 
and v(X) such that 



( u )t = ( v )\ 



(19) 
(20) 



The relationship between U and V (or between u and v) 
is highlighted in Appendix | A 1| Sections HI C and HI D| 



present the necessary Fourier analysis details. 

We will always avail ourselves of the simplication in 
equation (19). Accordingly, throughout the rest of the 



paper, we focus exclusively on torus averages over x & n d 
time averages over Mino time A with the understanding 
that they may merely be computation-friendly proxies for 
averages of different but related functions over different 
torus or time coordinates. 



Torus averaging and Fourier analysis of 
tor us- functions 



We will call a torus function f(x; £) any rule that as- 
signs a complex number to every point on a phase space 
torus. Note that £ specifies both the torus function and 
the phase-space torus that serves as its domain. Usually, 
we will be discussing properties of torus functions evalu- 
ated at some definite value of £. We thus omit the ex- 
plicit dependence on the orbital parameters £ for brevity, 
except where it might lead to confusion. 

We assume that every such torus-function is continu- 
ous and diffcrcntiablc in all its arguments (including £). 
We also require it to be single- valued on the torus, which 
implies 2n periodicity in each of the angle variables: 



f(Xr,Xe) = fiXr + 27r,xe) = f(Xr,Xe +2tt) 



(21) 



Like any function that is independently periodic in 
two independent variables, a torus-function has a dou- 
ble Fourier series representation 



f(x) = Y. A ^ ( 



-™Xr e -ikxe 



(22) 



k . n 



with the Afc n 's given as usual by 

Akn = — 2/ d Xr d Xe .f(Xr,Xe)e +m ^e+ tk *° . 
(2tt) Jo Jo 

(23) 

In order to distinguish them from another set of double- 
indexed quantities we introduce later, we will refer to the 
A^s as spatial Fourier coefficients or torus Fourier 
coefficients. 

We now define, in the usual way, the following useful 
quantities. The torus average of / is 



1 



27T rl-K 



(24) 



(fix)) ^=—~2 d Xr / dXfff(Xr,X0) 

[Zir) Jo Jo 
= Aio ■ 

The torus averaged Fourier power of / is 

V ^=7^2 I d Xr I dxe\f(xr,xe)\ 2 ■ (25) 
(2ir) Jo Jo 
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By Parseval's theorem, the torus-averaged power must 
also equal 



(26) 



k,n 



The 2D power spectrum of / is the contribution to the 
torus-averaged Fourier power from each pair of spatial 
frequencies or wavenumbers (n r ,ng). Note that because 
the period in each of the Xr and xe directions is 2ir, the 
corresponding fundamental spatial frequencies are K r = 
Kg = 1, so we see power only at integer lattice points 
(k, n) in the 2D wavenumber space. 

All statements above are standard results from the 
Fourier analysis of functions on a compact 2D spatial 
domain. They apply equally well on resonant and non- 
resonant tori. 



D. Time averaging and Fourier analysis of 
time-functions 



We can evaluate any torus function along a curve 



(12 1 on its associated torus that corresponds to an or- 
bit. Since each orbit is specified by its initial position 
Xo on the torus, each torus function naturally induces 
a 2-parameter family of time functions, one for each 
(Xr >X8 ) P au "- Time functions, then, are grouped into 
5-parameter families - 3 parameters £ to specify a torus, 
and 2 parameters xo to specify an orbit on that torus. 
We denote one memeber of such a family as 

/ (jcWlXo]£j ■ We will sometimes omit the explicit xo 

and ^-dependence of a time function and simply write 
/(A), again except where clarity would suffer. Through- 
out this paper, we adopt the notational convention that a 
time function and the torus function from which it is de- 
rived are denoted by the same symbol (/, in the examples 
so far). 

For time functions, non-resonant and resonant tori 
must be treated separately. 



1. Non-reosnant tori 

When Q()/Q, r is irrational, every k,n pair leads to a 
distinct frequency 



k,n G Z 



(27) 



Such a biperiodic time-function is not periodic: it is 
bounded on (—00,00), but there is no finite time inter- 
val over which the function exactly repeats itself. Still, 
every biperiodic function has a unique discrete Fourier 
representation [24] 



/(A;xo) 



(28) 



Note that the harmonics are not equally spaced in fre- 
quency. The temporal Fourier coefficients A n k-\ (which 
we have named suggestively) are given by the limit 



Akn-,\ = lim t 



b+A/2 



6-A/2 



d\f(X;xo)e l{nnr+kne)X 



(29) 



which exists and is independent of b |24j (henceforth, we 
set 6 = for convenience). Equivalently, we could say 
that the Fourier transform of /(A; xo) consists of a series 

7(fi) = Y, A kn^Hn-(nn r + kn e )) (30) 



of delta-function impulses unequally spaced in frequency. 

Paralleling the Fourier discussion of torus functions, 
we now define the time-averaged value, time-averaged 
Fourier power, and the ID power spectrum of a time 
function associated with a non-resonant orbit. Biperi- 
odic functions offer no single period over which to time- 
average in a natural way. Given the existence 11 of ex- 
pressions like (29), averaging over all time seems like a 



sensible choice. The theory of almost-periodic functions 
states that such an infinite time-average indeed exists 
[231. so we define 



(/0;xo))a = lim i 

A— ^00 A 



dA/(A;xo) 



(31) 



We will refer to (/), simply as the time- average of / 
rather than as the infinite time-average value, as it is 
sometimes called. Comparing equations (29) and (31), 
the time-average equals 



(/(A;xo)) ; 



An 



(32) 



We define the time-averaged Fourier power of /(A) as 
a special case of (31 1 by 



Vx = A lim T 

A— >oo A 



rfA|/(A;xo)r 



(33) 



Parseval's theorem also applies to biperiodic time- 
functions 25 , so 



fcn:A I 



(34) 



k :n 



11 Technically, the logical presentation of Fourier coefficients and 
time-averages for biperiodic (or more general multipcriodic) 
functions goes in the reverse order. First, the existence of the 
infinite time-average in | |29| l is established for a biperiodic func- 
tion /(A). The existence of the Fourier coefficients in ( |29| then 
follows from the existence of the average value and the fact that 
/(A)e* < - nf2 *" +fef28 ' A is also biperiodic. We have chosen this order 
to parallel the presentations in Sections IIII CI and IIII D 21 
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That fact allows us to define a ID power spectrum for 
/(A) as the contribution to the time- averaged power from 
each temporal frequency ft. The graph of |^4fcn ; A | 2 over 
the ID f2-space would show power only at the discrete 
and unequally spaced set of frequencies (27). 



The question now is how to evaluate these time av- 
erages in practice. Though equation (29) defines the 
Akn-.x's, such integrals over infinite intervals divided by 
infinite quantities do not lend themselves to simple eval- 
uation, either analytically or numerically 12 . 

To compute the temporal Fourier coefficients of 
/(A;xo)j we must instead proceed circuitously. Con- 
sider the spatial Fourier representation (22 ) of the torus- 
function fix) evaluated along the orbit (12), which yields 



/(A;Xr o ,X0 Q ) = ^A kr , 

k , n 



(35) 



By uniqueness 13 of the Fourier representation of /(A; Xo), 
and comparing equations ( 28 ) and ( 35 1 , we conclude that 
the temporal Fourier coefficients of /(A; Xo) an d the spa- 
tial Fourier coefficients of f{\) are related 14 by 



.4 



/»•//: A 



= A 



-inxr +kxe 



We note that each temporal coefficient differs from the 
corresponding spatial coefficient in ( 23 1 only by a com- 



plex phase determined by the initial conditions xo of the 
orbit. Consequently, their magnitudes are identical, re- 
gardless of the initial position of the orbit: 



\A kn ;X\ = \A kn \ , VxoeT 



EL Z Q 



(37) 



This fact is consistent with the ergodic property of these 
orbits. Every orbit eventually comes arbitrarily close to 
every point on the torus, so shifting initial conditions 
leads to a new orbit that is arbitrarily close to a time 



12 To evaluate equation j29j | numerically, larger and larger values of 
A would be required before converging to some accuracy. This is 
computationally impractical because such a process will in gen- 
eral converge extremely slowly. Thus, as the size of the integra- 
tion interval grows, so will the required number of evaluations of 
the integrand, a particularly problematic development if the in- 
tegrand is expensive to calculate. Moreover, the prefactor of 1/A 
can eventually become so small that there is loss of significance 
in the final answer, thus compromising accuracy. 

13 The set of complex exponential functions e _lf2A for all f! are a ba- 
sis for absolutely integrable functions on the space A £ (— oo, oo). 
/(A; xo) inherits absolute integrability from the associated torus 
function f(x), which has a spatial double Fourier series represen- 
tation and thus is absolutely integrable by assumption. Equation 
| |28| l is therefore a projection onto the complex exponential basis, 
and projections onto basis sets are unique. 

14 We denoted the fcnth temporal Fourier coefficient by A^n-x in 
anticipation of its close relationship to the fcnth spatial Fourier 
coefficient A^ n of the associated torus function and added the A 
subscript to remind us of when we are dealing with spatial vs. 
temporal Fourier coefficients. 



translation of the original orbit. And, of course, time 
translations only affect the complex phase of temporal 
Fourier coefficients. 

If we know the torus function fix), its spatial Fourier 
coefficients A kn can be computed by any number of ef- 
ficient numerical routines, without any of the difficulties 
that beset computation of the A kn .\S via direct evalu- 
ation of the definition (28). This fact, combined with 
equation (|36| 



leads to t 



re only practical recipe for com- 
puting the Akn-,\S of the orbit with initial position X0j 
namely to compute instead the A^'s and then use equa- 
tion (36). Rcf. [55] introduced just such a technique in 



the specific context of functions of Kerr geodesies. 

All the other quantities mentioned in this section are 
likewise determined from their torus function counter- 
parts. From equations (32) and (36), ^4qo ; a = ^oo- We 



thus conclude that on a non-resonant torus, the time av- 
erage of /(A) equals the torus average of its associated 
torus-function fix)- Moreover, since this is true for ev- 
ery time function on that torus, the time average of such 
a function is independent of the initial condition xV 



(/(A;xo)) A 



(38) 



(36) Likewise, equations (37) and (26) imply that, on a 



non-resonant torus, the time-averaged Fourier power of 
/(A; xo) equals the torus-averaged Fourier power of fix) 
for every xo'- 



(39) 



Equations (26), (37) and (39) together imply equation 
(34) 15 . By extension, the ID power spectrum of /(A) 



can be derived from the 2D power spectrum of fix) by 
mapping wavenumber pairs to frequencies using eq. ( |27[ ) . 

All of the above relationships between torus-function 
quantities and those of any biperiodic time function in- 
duced via ( 35 ) are well-established and well-known in the 
literature on almost-periodic functions [24, 25 and on in- 
tegrable Hamiltonian systems [TS]. Many of these facts, 
however, are used but not so clearly delineated in this 
way in the literature relating to EMRI calculations. We 
have gone through the trouble of including them here not 
only for completeness and clarity but also to emphasize 
that we can only execute the above recipes if we know 
the corresponding torus function fix)- 

This leads us to a crucial observation. If all we know 
is /(A), either as some closed- form expression in terms 
of A or as a numerical time-series, there is no practical 
scheme for computing its temporal Fourier coefficients 
Akn-x, even though those coefficients are perfectly well- 
defined. In addition to the initial conditions xo associ- 
ated with /(A), we must also know the torus-function 
fix) ( or & f least its spatial Fourier coefficients Af. n ) in 



This proves Parseval's theorem for biperiodic functions. 
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order to compute the Ak n ;\s. We summarize the impli- 



cations of this fact for flux balancing in section III E 



All of the equivalences noted between torus-function 
quantities and their time-function counterparts followed 
from the assumption of non-resonance. On resonant tori, 
all of these equivalences break down, as we now show. 



2. Resonant tori 

Unlike time functions evaluated along non-resonant or- 
bits, time functions on resonant orbits are singly periodic, 
with (possibly very long) period Ap and corresponding 
fundamental frequency ilp = 2n/Ap. The single peri- 
odicity of time functions of resonant orbits means that 
all frequency-domain quantities have straightforward and 
familiar definitions. 

Any time function evaluated on a resonant orbit has a 
Fourier series representation 



ij^lp A 



whose coefficients 16 are single-index objects 

1 /"Ap/2 



C 



A, 



P J-hp/2 



dXf(X;xo)e + ^ npX 



(40) 



(41) 



Like the Ak n .\S, each Cj- t \ varies with the initial con- 
dition xq. Unlike the Ak n -,\s, the Fourier transform of 
/(A; xo) is a sequence 



(42) 



of equally spaced delta-function impulses in frequency 
space. 

In the resonant case, we define the time average of /(A) 
straightforwardly as 



(f(X;xo))x 



l 

a7 



dXf(X;xo) = Co ; a 



Likewise, we define the time-averaged power as 



^a(xo) 



A, 



p j p 



dX |/(A;xo)r 



(43) 



(44) 



By Parseval's theorem, the time-averaged power is also 
given by 



(45) 



Periodicity of /(A) implies that the integral in equation 

has the same value taken over any interval of length Ap. We 

choose the symmetric interval [—Ap/2, Ap/2] solely for aesthetic 



To flush out how the Cj-\s relate to the spatial A^n's 
and to the initial condition X0i we begin, as in the non- 
resonant case, by inducing a time function ( 35 1 from a 
torus function. In the resonant case, the frequency ra- 
tio Qg/Q r is a rational number p/z, where p and z are 
relatively prime and p > z. (In terms of integers in the 
definition q r g = w+~,p = (w + l)z+v.) The individual r 
and 9 frequencies and periods are therefore related to the 
fundamental frequency and total period of the periodic 
orbit by 



O r — z£2 p 

fig = pflp 



and 



A r = 



An 



Ap 

z 

Ap 
P 



(46a) 
(46b) 



(47a) 
(47b) 



As a result, all kn combinations for which 

nz + kp — j (48) 

lead to identical frequencies 

ntt r + kQg = nzflp + kpflp = jQp (49) 

in the arguments of the exponential functions on the RHS 
of equation ( 35 1 . 



The selection rule (48) maps every kn pair to some j. 



By the uniqueness of Fourier represenations, we conclude 
from equations ([35l and (|40l that 



Cr,x(xo)= £ A 



-i(riXr +kXe ) 



(50) 



k,n: 
nz+kp—j 



Note that equation ( 50 ) is really only a summation over a 



single index since the value of k in any term is determined 
by the value of n and the (fixed) value of j. 

It is tempting to rewrite each term on the RHS of 
equation (50) as Ak n -x, mimicking the notation for the 
non-resonant temporal Fourier coefficients. We refrain 
from doing so because we seek a clear distinction between 
spatial and temporal Fourier coefficients, and temporal 
double-index coefficients are not defined in the resonant 
case [251 HZ]- Fourier representations are unique, so the 
familiar single-index representation ( 40 ) is the only such 



projection of /(A) onto a set of mutually orthogonal basis 
functions. If we were to write an expression like (28 1 on a 



resonant orbit, the different harmonics on the RHS would 
not all be orthogonal, and we would not have a bona fide 
Fourier expansion in hand until we collapsed all terms 
corresponding to the same frequency into a single term. 
Equation (36) implied that, on non-resonant orbits, 



several quantities one can compute for a time-function 
f(X;xo) turn out to be independent of xo : the magni- 
tudes of its Fourier coefficients, its time-averaged value, 
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its time-averaged Fourier power, and its power spectrum. 
In contrast, equation (50 1 implies that, on resonant or- 



bits, each of those quantities does depend on the initial 
condition xo- Each Cj-\ is a sum of spatial A^nS with 
Xo-depcndcnt phases rather than just one such term (cf. 
equation (36)), so both the magnitudes 



\C j; x\ (Xo) 



-i{nxr +kxe ) 



k,n: 
nz+kp—j 



and time-averaged value 

(/(A;xo)) a = 0> ; a(xo) 

= ^2 Akne~ i( - nXro+kxea) 



k,n: 
nz-\-kp—0 



A 



oo 



E Akn< 



-i(nx rQ +kxe ) 



n^0,M0: 
nz+kp— 



(51) 



(52) 



retain xVdependence. The squared magnitudes 



\CAXr ,X6 )? = £ Ak n e-^ +k X(>0 ) £ Al, n A n '^o+«xe ) 

kn: k'n': 
j—nz-\-kp j=n / z-\-k'p 

= J2 \ A kn\ 2 + £ A kn Al, n ,e-< n - n ')^e-< k - k ')^ 



(53) 



k=k ,n=n : 
j=nz-\-kp 



k^k ,n^n : 
j—nz+kp, 
j—nz+k p 



also depend on \o through cross terms, and the time- 
averaged power and power spectra inherit this depen- 
dence via ( 45 ) . Figure [6] illustrates this point for the test 
function r cos 6. 



Via its xVdependence, equation (52 1 defines a torus 



function in the variables Xr , Xo ■ Complex exponentials 
have a zero average value, so averaging that torus func- 
tion over all x rn , xb kills every term in the summation on 
the RHS of (52), leaving only A no . But A 00 is the torus 



averaged value of the associated torus function f(x)- We 
conclude that the torus-average over all initial condi- 
tions of the time average of a time function equals the 
torus-average of the underlying torus-function. An iden- 
tical argument applies if we torus average the squared- 
magnitudes (53) of the coefficients over all xo and, by 



extension, if we likewise torus-average the time-averaged 
power (44). 



The upshot is that the parallels between torus func- 
tions and time functions obtained in the non-resonant 
case break down in the resonant case: 



(/(A;*o)) A ?M/(Xr,X9)>* 



(54) 
(55) 



However, torus averages over initial conditions and torus 
averages of time averaged are equal for both a function 
/ and its Fourier power: 



(Vx(xo)) 9 = V x . 



(56) 
(57) 



To clarify, equations ( 54 ) and ( 55 ) state that the time 



average of a time function and the torus average of its as- 
sociated torus-function are not identically equal as they 
are in the non-resonant case. That does not, of course, 
preclude the possibility that the two could be circumstan- 
tially equal for some particular choice of initial condition 
Xo- In fact, for real- valued functions /, the mean- value 
theorem guarantees that (f{\;xo))\ = (/(x)}-? f° r at 
least one x|J wt € T|.. In general, / will be complex- 
valued, and we have no such guarantee. The time- 
averaged Fourier power, however, is strictly real, so there 
is at least one x^ vt such that V\{xT*) = ^ We ex- 
plore some implications of this fact for adiabatic EMRI 
calculations in Section IY1 
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FIG. 6: Top: The Fourier power spectrum of the function 
r cos 9 for a q r e — | periodic orbit for two different sets of ini- 
tial coordinates, (r = r a , 9o = f ) and (r = r a ,0 = m4n ), 
but the same sets of orbital parameters, a = 0.9, E = 
0.954788, L z = 2.65115 and Q = 0.944969. Bottom: The 
magnitudes of some spatial Fourier coefficients for the same 
orbits in the top panel. 



E. And the winner is. . . torus averaging 

To summarize, time averaging is equivalent to torus av- 
eraging for non-resonant orbits. Furthermore, torus av- 
eraging is the only practical recipe for computing Fourier 
coefficients and so torus averaging is the explicit compu- 
tation instituted in practice. 

However, time averaging is inequivalent to torus av- 
eraging for resonant orbits. Thus torus averaging wins 
out for two reasons. First, torus averaging along a 
resonant torus crucially washes away any xo positional 
dependence, while time averaging does not. The Xo~ 
dependence violates the spirit of averaging, namely to 
remove all dependence on the fast variables. Second, and 
more seriously, though the time-averaged equations are 
continuous in xo f° r fixed £, they are in general discon- 



tinuous in £ for fixed xo- The situation will resemble 
that of Thomae's modified Dirichlet function 

if x is irrational 
Dm(x) = \ \ if x = p/z, with p and z coprime , 
if x = 

(58) 

which is continuous at the irrationals, discontinuous on 
the rationals, and nowhere differentiable 17 [29]. Prag- 
matically speaking, even if a set of ODEs with such 
pathologically discontinuous and non-differentiable equa- 
tions had a solution, it is unclear how one would numer- 
ically integrate them. Furthermore, the continuity fur- 
nished by torus-averaged fluxes is absolutely essential for 
the proper construction of a grid through which adia- 
batic trajectories are to be interpolated, as discussed in 
the introduction. 

In short, torus-averaged equations are well-behaved, 
while time-averaged equations lose the continuity and dif- 
ferentiability that guarantee the resulting equations are 
well-posed and have unique solutions, the very basis of 
every standard numerical integration scheme. 

The arguments made in favor of torus averaging ap- 
ply to the radiative approximation, based on an aver- 
age of the dissipative piece of the local self-force on the 
inspiraling particle. But when it comes to flux-balance 
as a specific implementation of the radiative approxima- 
tion, this now leaves a logical gap. As acknowledged in 
[TJ [11] , the flux-balance arguments that allow the nonlo- 
cal fluxes of conserved quantities to be used as proxies 
for the local dissipative self-force have been derived on 
a time-averaged basis and under the assumption of non- 
resonance. Since time and torus averages agree for non- 
resonant orbits/tori, the time-averaged nonlocal fluxes 
are still good proxies for the torus-averaged local dissi- 
pative self-force in the non-resonant case. 

What has never been made explicit is whether flux- 
balance is also valid in the resonant case, on either a 
time averaged or torus-averaged basis. We resolve this 
issue now: time averaged flux-balancing may not be true 
on resonant orbits in general, but that will be irrelevant 
since it will be true on a torus-averaged basis. Mino 
showed that, under the assumption of non-resonance, the 
time-averaged fluxes of £ at infinity and the horizon fur- 
nish proxies for the time averaged RHSs of equations 
(18). But then, by the arguments we have heavily ex- 



ploited, the corresponding torus-averaged versions must 



17 Such a function certainly seems unphysical. It violates the hy- 
potheses of continuity and differentiability in all arguments re- 
quired by the theorems bounding the error in a solution to time- 
averaged equations with almost-periodic dependence on time 
|17| . More pathologically, it violates the hypotheses for the well- 
posedness of an initial value problem and for the existence of 
solutions to systems of ODEs 28]- Such a function would not be 
Riemann integrable, and it is unclear whether a system of such 
functions would even be Lebesgue integrable. 
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also be equal. And since torus averages are insensitive to 
the resonance or non-resonance of the underlying torus 
and are continuous in £, the flux-balance prescription 
is valid in a torus-averaged sense for all orbits. If the 
torus-averaged fluxes were not good proxies for the torus- 
averaged local equations only at resonances, a discrete set 
of measure zero, then they could not be continuous in £. 
But we have shown torus averages are everywhere con- 
tinuous. De facto, then, Mino's argument establishes the 
validity of torus-averaged flux-balancing generally. 

It would thus seem that both flux-balancing as a 
general procedure and its specific implementation in a 
frequency-domain application of the Teukolsky formalism 
treat non-resonant and resonant tori equally, as stated in 
[7]. Flux-balancing is, in fact, thusly impartial, but in- 
terestingly, the Teukolsky formalism is not. As we will 
show below, a TB torus-averaged flux calculation can 
achieve computational savings of an order of magnitude 
or more on low-order resonant tori that are simply not 
available on non-resonant tori. These efficiencies follow 
from a simple observation about the Fourier integrals a 
TB code must evaluate and are independent of the spe- 
cific implementation in code of the Teukolsky formalism. 
Thus, rather than disfavoring resonances, as has been 
commonly assumed, the Teukolsky formalism actually 
shows favoritism for resonances, and properly leveraged, 
that favoritism can substantially accelerate adiabatic in- 
spiral calculations. 



IV. COMPUTATIONAL SAVINGS ALONG 
RESONANCES 

We now explain computational efficiences that exploit 
resonant tori. Although specific to the Teukolsky formal- 
ism, the computational expedience can be understood 
without all details of that formalism. We simply assert 
some features and formulae from a TB flux calculation 
that we require to make our argument. For reference, 
Appendix [B] gives a somewhat more detailed overview of 
the Teukolsky formalism and offers at least skeletal devi- 
ations of the formulae listed below. For a fuller treatment 
of the Teukolsky formalism, which is beyond the scope of 
this work, we direct the reader to the references listed in 
Appendix [B] 

As a reminder, our argument is specific to frequency- 
domain Teukolsky calculations and corresponding codes. 
In the context of the EMRI problem, such codes com- 
pute a combined multipole and Fourier decomposition of 
the metric perturbations at infinity and the black hole 
horizon due to a geodesic source. 



on non-resonant orbits 18 . By the arguments of Section 



HID 1 and Appendix [AJ this is equivalent to the torus- 
averages of the fluxes on all tori over the torus coordi- 
nates 7. We therefore report those same expressions here 
as torus-averaged fluxes. For E and L z , those expressions 
are [3 [14] 



dE\ 
~dt L 

I 7 



n 



00/H 
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H/c 



Imkn 

Imkn ^^mkn 

-j/H 
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J lmkn 



(59a) 
(59b) 



Based on Mino's argument in [TTJ, Refs. [THIGH] worked 
out the corresponding expression for the time-averaged 
Q flux for non-resonant orbits, which we also report as 
the torus-averaged flux 



— ) = -2 (a 2 E cos 2 6) J — \ 



dt / 7 - x ^\dt 

H/oo 



(59c) 



Imkn 



7 oo/H 
J lmkn 



The prefactors in the first two lines of (59c 
puted only once for the entire torus 
and 



are com- 



(59b) 



Thus, substituting 
into the RHS of fl59cl and 



equations ( 59a 
combining like terms with those in the summation of the 
last line, the flux for Q has the same general form as the 
fluxes of E and L z have. Our savings arguments will be 
based on that form, so although we will speak about E 
and L z for concreteness, those arguments will apply to Q 
as well. Appendices |B] and [C] summarize the derivations 
of these expressions. 

Before proceeding with those arguments, we clarify the 
notation in equations (59). First, the apparent discrep- 



ancy between the ordering of the H/oo superscripts on 
the left- and righthand sides of the equations is not a ty- 
pographical error. On the LHS, the superscript denotes 
fluxes at the black hole horizon and radial infinity, re- 
spectively. The somewhat backward notational choice to 
have the fluxes at infinity depend on a quantity labeled 
^irnkn an d the horizon fluxes on Zff nkn is, at this point, 
ingrained in the literature. To maintain a modicum of 
notational uniformity, we have labeled the weighting fac- 
tors Oi H l°° with the same backward superscript conven- 
tion. The exact form of those weighting factors will not 
concern us. What matters for our purposes is that every 



factor 



l lmkn 



for the fluxes at infinity is equal to 1 and 



that every factor a^ fcn for the fluxes at the horizon is 



A. The fluxes of E,L Z ,Q 



The fluxes of the conserved quantities £ are usually 
reported as quantities averaged over coordinate time t 



18 See Appendix [c] for the time-averaged fluxes from resonant or- 
bits. 
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real and depends on k,n only through w/- n . All the ar- 
guments to follow apply equally to fluxes at infinity and 
at the horizon. We borrow the notation * from Ref. [2] 
to denote either of H/oo. 

Continuing, the indices I, m are standard multipole in- 
dices 19 , with I > 2, — I < m < I. Our argument will focus 
on the Fourier analysis of each I, m term individually, so 
that, unless explicitly stated otherwise, l,m are taken to 
be fixed everywhere in this section, while k, n each run 
from — oo to oo. The frequencies 



kua 



(60) 



are the combined harmonics ojkn of the r and 8 funda- 
mental frequencies (the coordinate time version of equa- 
tion (27)) and the fundamental azimuthal frequency lj v . 
Note that the integer m is both a multipole index and 
the relative contribution of lo v to each frequency uj mkn . 
Other than attaching itself as a label to frequencies in 
this way, however, m will not appear as a Fourier index 
in any sense below. 

Finally, Appendix [C] explains why we have written the 
fluxes as average values over the 7 torus coordinates men- 
tioned in Section [ill B| and in Appendix [A] We note here 
simply that if we seek adiabatic solutions in the form 
£ (t) (as opposed to £(A)), then the angle brackets in 
(18) should also be averages over 7, so that (59) have 



the correct form to be proxies for F{£). The representa- 
tions of the LHSs of the flux equations as averages over 
7 is otherwise irrelevant, since in light of (19), we will 



always seek equivalent and easier to compute x-averaged 
quantities. 



as Fourier coefficients of a torus function 



Equation (61) parallels the form of equation (23) from 
Section HI C but there is one critical difference. For 
fixed l,m, the function fi m . u {Xr,Xe) further depends on 
a continuous parameter uj that must be set to uimkn when 
evaluating Zf kn for a given multipole mode. Postponing 
for the moment any details of the function f* m . u (Xr) Xe) 
or its derivation, we remark that this dependence on the 
coordinate time harmonic frequencies of the source as an 
external parameter persists despite the fact that equation 



(61 1 is a spatial Fourier integral. 

Thus, for fixed I, m, the Z* mkn are not the Fourier co- 
efficients of a single function but rather isolated Fourier 
coefficients of several different functions 20 . On a non- 
resonant torus, every k, n pair leads to a different value of 
Umkn , and every coefficient computed has a distinct func- 
tion f* m . u (Xri Xe) m the integrand. On a resonant torus 
with associated frequency ratio LOg/uj r = Q,g /Cl r = p/z, 
all k,n pairs that satisfy the selection rule (48) for the 
same j lead to identical values of fi m fcnj and some coef- 
ficients with different values of k, n will share the same 
integrand function f^.^ixr, Xe)- The practical implica- 
tions of this asymmetry for a TB flux calculation consti- 
tute the basis of our savings argument. 

In anticipation of later arguments, we also note that 
for each fixed value of the pair I, m, the resulting doubly 
infinite sum over k,n in (59) has the appearance of a 



torus-averaged Fourier power in the sense of Section III C 
with the identification 



\A kn \ 2 = prefactor x \Z^ mkn \ 2 



A kn — \J prefactor x Z* r 



(62) 



rnkn 



The prefactors in front of \Zi mkn \ turn out to be real- 
valued and non-negative for all values of the indices, so it 
is valid to subsume them into some new coefficients A kn . 



With these preliminaries out of the way, we are ready 
to list the features of the RHSs of ( 59 ) that we will need 



for our savings arguments both in this section and in 
Section [V] For our principal argument, what matters is 
that for fixed I, m values, each Z* mkn takes the form of a 
Fourier coefficient of some torus function, 



7* = 

Imkn 

" d Xr d X e e^ +k ^fL- M =^ mkn (Xr,Xe) 



(61) 



This form of the Z* mkn y s associated with a geodesic 
source of arbitrary eccentricity and inclination is detailed 
in several references (see, for instance, [3J Q31 HS1 HH [2E1 
130] ) and summarized in Appendix [B| 



C. Recycling computations between Fourier modes 

The complex- valued quantities Zt k are the backbone 
of a frequency-domain Teukolsky calculation, and a code 
that implements such a calculation spends by far the 
lion's share of its CPU budget on computing them. To 
explain how resonances can be leveraged to optimize that 
budget, we must look a bit more closely at the integrand 
functions /,* m . w (xr,Xe)- 

The main ingredients in f* m . u (Xr>X9) are two sepa- 
rate functions that have the same sort of uj dependence 
described above: a radial Teukolsky function R* m . u (r) 
and a spin- weighted spheroidal harmonic -2&f" (0). We 
imagine re-expressing the former as a torus function of 
Xr alone and the latter as a torus function of xe alone 



The values I = 0, 1 are not relevant in GW calculations, for which 
the lowest non-vanishing moment is the 1 = 2 quadrupole. 



This is part of the reason why we cannot compute all the Z£ mkn 
coefficients for a given I, m at once with, for instance, a 2- 
dimensional Fast Fourier Transform (FFT). 
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but will continue to write them as functions of r and 9, as 
they are in the rest of the literature. fi m - u {Xr , X&) con- 
sists of a somewhat messy assortment of terms and fac- 
tors involving these two functions, several of their deriva- 
tives, the coordinates and velocities of the particle (both 
of these are absorbed into the torus coordinates Xr , Xe ) , 
and other elementary functions. 

Each of Rf m . u and - 2 S?£ (9) satisfies an ODE that de- 
pends on /,m and w in a nontrivial and partly implicit 
way (see Appendix [b]). No simple closed-form solutions 
to these equations exist that make the functional depen- 
dence of the solutions on those parameters explicit. As 
a result, for every distinct set of values (I, m, u), those 
ODEs must be solved from scratch to obtain the numeri- 
cal representations of R~i rn . M and -zSff^ needed to evalu- 
ate the integrand. Particularly in the case of -R* OT . W , this 
operation is computationally costly. 

Schematically, then, one calculates each Zf mk for fixed 
I, m via the following steps: 

1. Determine the frequency lo = w mkn 



2. Obtain a representation of -2S' 



I in 



3. Obtain a representation of R* m . u for * = H, 00 (this 
step requires first determining an eigenvalue of the 
-aS^ODE) 

4. Evaluate f* m - u (Xi) at whatever abscissae Xi are re- 
quired by the specific numerical integration algo- 
rithm chosen 



5. 



Compute whatever weights to, the integration al- 
gorithm may require for those function values and 
tabulate the integral ([6ll as £\ iVifi m . u (Xi) ■ 



On a non-resonant torus, each k, n pair produces a dif- 
ferent answer to step 1 and requires the execution from 
scratch of all the remaining steps as well. On a reso- 
nant torus, in contrast, the k, n pairs can be grouped by 
a common value of j in the selection rule ( 48 ) . Steps 



1-3 need only be performed once for an entire j-group. 
Depending on the integration algorithm selected, steps 
4 and 5 may also only need to be performed once or a 
small number of times per j-group, with a total number 
of reusable function evaluations set by the Zf mkn in the 
group requiring the greatest number of sample points to 
attain some target accuracy. We will make the reason- 
able assumption that steps 1 and 5, even if done several 
times per j-group, are a small fraction of the total cost of 
evaluating all the coefficients in that group, and we will 
take the cost of steps 2-4 as an estimate of the total cost 
of computing any single coefficient. 

Consider now evaluating all the Zf mkn on a low-order 
resonant torus with u)e/uj r — Qg/£l r = p/z = 1 + q r g and 
on a neighboring non-resonant torus with nearly identical 
orbital parameters. By the continuity of the Z^ mkn with 

respect to £, coefficient values will be nearly identical on 
those two tori. The integer values of n max , fc max deter- 
mined should also be identical or nearly identical on the 



two tori (we assume for simplicity that they are identi- 
cal). Let A/2 and A/i denote, respectively, the number of 
separate times steps 2-4 above must be executed on the 
non-resonant torus and resonant torus. To make a more 
apples to apples comparison, one can instead let A/2 rep- 
resent the total number of distinct executions of steps 2-4 
on the resonant torus if the resonance of that torus is not 
acknowledged from the outset. Roughly speaking, gen- 
erating all the Zf mkn with |n| < n max , |fe| < fcmax on the 
non-resonant torus will require A/2/A/1 times more com- 
putation than it will on the neighboring resonant torus. 
Symmetries in the underlying equations imply that the 
value of Z*/ m -<, k <.,_ n \ is uniquely determined by the 
value of Z* mnk . Thus, in practice, one of the indices n 
and k can be restricted to run over only nonnegative val- 
ues, and the value of A/2 /A/i should take that fact into 
account. 

Figure [7] estimates the savings factor A/2 /A/i for reso- 
nant tori with various values of q r $ and for several rep- 
resentative hypothetical values of n max and fc max consis- 
tent with the reported performance of the TB code for 
arbitrary eccentricities and inclinations described in Ref. 
|14) 21 . For simplicity, we have taken n max = fc max . 

We can see the following trends in the histograms. 
First, for a given fc max , n max , if we fix the value of p and 
increase z or vice versa, the savings factor drops. Thus, 
the savings factor is largest when both p and z are as 
low as possible. The greatest savings (over an order of 
magnitude) accrue when z = 1. Second, the larger the 
values of 

^ma^i^maxi i-e. the more slowly converging the 
expressions for the fluxes, the greater the savings fac- 
tor for a given p/z. Generally speaking, the most slowly 
converging fluxes are for orbits with moderate to high 
eccentrities [2J [TU] , which typically have higher associ- 
ated values of q r e since they are closer to the separatrix 
between plunging and non-plunging motion [S]. Thus, 
for instance, a rough approximation of the true savings 
factor in the top two panels of Figure [7] would be given 
by a roughly horizontal or slightly downward sloping line 
connecting the histogram bar with the lowest q r g at the 
lowest n ma x with the highest q r g at the highest n max . 
A good rough predictor for the expected savings would 
thus be the z value of a torus, yielding a savings factor 
of ~ 30 for z = 1 and ~ 7 for z — 6. The lowest sav- 
ings factor on that graph of ~ 3 (corresponding to the 
not-so-low-order resonance with p/z — 25/6) is nothing 
to sneeze at, and more typically the savings factor from 
acknowledging resonance would appear to be around an 
order of magnitude on average. 

While a detailed audit of comparative cost would have 



21 We estimate n max and fc ma x based on the code in 1141 rather 
than the similar code in 1231 only because the truncation rules 
used in [14) are more amenable to direct cost comparison with 
our proposal. Both codes seem to need to compute a total num- 
ber of modes of similar order of magnitude to achieve high flux 
accuracy, and both are apt to profit from our proposal. 
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FIG. 7: The three histograms show the average number 
of Z* mnk coefficients that pertain to a single frequency on 
a resonant torus, a number that corresponds to the sav- 
ings factor JV2/A/1. We show the savings factor for a vari- 
ety of q r g geodesies and a variety of rc max and fc max . Top: 
q r e = Integers. Middle: q r g is a variety of values all with 
the same denominator, z = 6. Bottom: q r g is a variety of 
non-integer values all with p = 11 but different z. 



to be done on a code-specific basis, the potential payoff of 
these observations makes a case for testing our proposal 
in existing codes. 



D. Numerical EMRI grids 

Even if we stipulate that fluxes can be computed 
more efficiently on low-order resonant tori than on non- 
resonant tori, is this fact necessarily useful? After all, 
points in £-space corresponding to resonant tori, let alone 
low-order ones, are already a measure zero set, so to con- 
struct the inspiral curve £(t), wouldn't the RHSs of the 
adiabatic ODEs have to be evaluated in general (and, 
formally, infinitely more often) on non-resonant tori than 
on resonant ones? Interestingly, while the answer to that 
question is "yes" , the actual calculation of TB fluxes itself 
need only ever be done on resonant tori, and at least pre- 
dominantly (and possibly exclusively) on low-order ones, 
at least for the foreseeable future. 

The reason has to do with the absolute computational 
cost of those fluxes, even on resonant tori. Simply insert- 
ing a TB frequency-domain flux routine into the RHS 
of, say, a standard Runge-Kutta ODE solver to gener- 
ate inspiral curves in real-time is untenable, even with 
a large number of processors at one's disposal to paral- 
lelize the TB calculation. Instead, solution of the adia- 
batic ODEs will proceed as follows. For each value of the 
black hole parameters, one would build a numerical grid 
of flux values on some dense mesh of points in £-space 
and then interpolate off of that grid to obtain the fluxes 
for arbitrary values of £. Once handed such a grid, those 
interpolated flux values would go into a standard ODE 
solver which could presumably generate inspiral curves 
very efficiently. The main expense to consider, then, is 
the construction of the grid. 

Our proposal is that such a grid should be built using 
exclusively resonant grid points. More specifically, we 
propose a hierarchical population of such a grid, begin- 
ning with the low-order resonant points and then increas- 
ing the order of the resonance (or just increasing z, if our 
loose conjecture about horizontal lines in the top pan- 
els of Figure [7] proves to be correct) until some requisite 
grid density is obtained to minimize interpolation error. 
Those grid density requirements may force the evaluation 
of fluxes on some higher-order resonances, but no reso- 
nant grid point (whether low- or high-order) will ever be 
more expensive to populate with Fourier flux data than 
a nearby non-resonant grid point will be. The worst-case 
scenario near certain locations in the space would be to 
break even by using a resonant versus a non-resonant 
grid point. Our hierarchical approach would seem to at 
least lower if not minimize the total computational cost 
of such a grid. 

We remark on two features of our proposal. First, 
the savings factor discussed depends only on the deci- 
sion to use resonant tori for TB calculations and no 
other implementation-specific features of that calcula- 
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tion. Thus, those savings will multiply any additional 
savings that may stem from other algorithmic improve- 
ments in any such implementation or from the availability 
of more or faster processors to perform the TB calcula- 
tions. 

Second, its efficacy has nothing to do with interesting 
physical effects that may occur in the neighborhood of 
resonant orbital parameter values during a real inspiral 
[2]. The flux-balance method, and in fact the adiabatic 
approximation in general, may fail to capture these ef- 
fects. Any such failure is immaterial to our argument, 
which rests not on physical properties of resonant tori but 
rather mathematical ones they have in specific relation to 
frequency-domain TB calculations. In other words, de- 
spite the fact that the adiabatic approximation might be 
least faithful to reality in and around resonances, leverag- 
ing resonances is nonetheless the most efficient means of 
attaining an adiabatic approximation for those regimes 
where it is likely to be faithful. 



E. Gravitational waveform snapshots 



As already argued, we are free to interpret the coef- 
ficients Z* mkn either as spatial Fourier coefficients of a 
torus function of x or °f a different torus function of 
7. As shown in Appendix |Bj the ^-function versions of 
the Z* mkn coefficients are used to build the Weyl scalar 
■04 at radial infinity, from which the two polarizations of 
the waveform h are constructed. These waveform "snap- 
shots" from geodesic sources [H| are useful for exploring 
how a known orbital motion impacts GW signals, and 
though they will quickly go out of phase with a true in- 
spiral signal, they are still likely to play a pivotal role in 
hierarchical searches for GWs from EMRIs. 

More specifically, with the 7-coefficients Zf mkn in 
hand, h can be reconstructed from the associated t- 



coefficients. By analogy to equation (36), we get 



7* = 7* f 

lmkn:t — lmkn c 



-i{n~f ro +k~(e ) 



(63) 



in the nonresonant case, and by analogy to equation (50) 
get 



E 

nz-\-kp—j 



-«(n7r +/E7e ) 



(64) 



in the resonant case if we know the initial conditions. 
Since the waveforms (or, rather, their Fourier represen- 
tations) depend on the Zf mkn coefficients, then like the 
fluxes, they will also probably need to be interpolated 
from a grid that stores the Z* mkn values themselves in- 
stead of or in addition to the fluxes. The same arguments 
made above for the fluxes thus cross-apply to waveform 
snapshots. 



V. SPECULATIONS ON FURTHER SAVINGS 

In this section, we sketch a speculative but tantaliz- 
ing possibility for further efficiencies in adiabatic EMM 
grid construction beyond those discussed in Section |IV| 
The idea centers around calculating time-averaged rather 
than torus-averaged fluxes on resonant tori. At first 
glance, that suggestion seems to fly in the face of ear- 
lier arguments that the RHSs of the adiabatic equations 
should be torus-averaged fluxes and that torus averages 
and time averages are not identical on resonant tori. The 
apparent incongruity disappears, however, in light of two 
facts: 

1. On any resonant torus, the mean value theorem 
guarantees that torus-averaged fluxes equal time- 
averaged fluxes on certain special orbits. 

2. For low-order resonances, those time-averaged 
fluxes are more accurate and cheaper to compute. 

The additional savings are beyond the cost benefit of 
incorporating the proposal of Section |TV} 

We substantiate these claims below in turn. We cau- 
tion the reader that, in contrast to the savings of Section 
|IV[ those discussed in this section may prove more elusive 
in practice because determining the special orbits men- 
tioned in step 1 above could prove so difficult as not to be 
net-beneficial. We discuss such limitations and suggest 
fruitful avenues of numerical investigation to help further 
reduce the cost of generating adiabatic inspirals. 



A. Using time-averages to compute torus-averages 

The time-averaged fluxes from a single resonant orbit 
do not appear elsewhere in the literature. As we explain 
in Appendix [Cj the arguments of Section |III D 2| imply 
that those fluxes are (note, these are single-index objects 
in j) 



dE\* 
dt/ t 

dt 



E a lmj 
Attuj 2 ■ 

Imj 



a lmj m 
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(65) 
(66) 



It remains to be shown whether the following would 

translate to ■ We restrict attention in this section 

to E and L z fluxes for the sake of exposition. 

As before, we assume fixed l,m in everything below. 
In the fluxes, the frequencies 



- ]U)p 



(67) 



the real- valued weight factors <x* m j, and the temporal 
Fourier coefficients 



7* 



1 

A7 



dXe^ fL (x(A;xo)) (68) 
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all become single-index quantities by the arguments of 
Section III D 2 We recall from equation ( 62 ) that the 



torus-averaged fluxes have the form of a torus-averaged 
power of some unspecified torus-function. Likewise, for 
fixed values of I, to, the time-averaged fluxes ( 65 ) and ( 66 ) 



have the appearance of a time-averaged Fourier power in 
the sense of Section HHP 21 with the identification 



|Cj;a| 2 =prefactor x \Zf mj . x \ 
C r . x = vVefactor X Z* mj . x 



(69) 



Each time-averaged flux, like any time-averaged 
Fourier power, is real-valued. Therefore, by equations 
(691 and (57 1 and the mean- value argument made at 

there exist initial positions 



the end of Section |HID2 

-*mvt,+ -*mvt,* 
Xq-E 'Xo ; L 2 



on the torus such that 



dE\ /^mvt,* 



~dT/ t ^ ) 



= / dE Y 



dt /- 



dL z 
lit 



(70) 
(71) 



Actually, there must be at least two continuous 1- 
parameter families of special initial values X(>e '* ( one f° r 
each of * = H , oo) and two such families for Xckl '* : an y 
two initial conditions that lie on the same orbit simply 
time-translate that orbit, and time-translation does not 
change time-averaged function values or time-averaged 
powers. 

None of the values X*cke'* and X*q-l '* nee d agree. Thus, 
if we sought to determine the torus-averaged fluxes in- 
directly by instead evaluating time-averaged fluxes, we 
might need to evaluate each coefficient Zf m j. x as many 
as four 22 times, once for each of the initial conditions 

X ,e and x ; Lz ■ 

We can, however, also apply the mean-value argu- 



ment individually to each real-valued 



7* 



In this 



case, we would obtain a sequence of special initial condi- 

2 



tions X^o-T'* that cause each 



7* 



to attain its torus- 



averaged value over all possible initial conditions. The 
different X(w '* would not necessarily agree for different 



values of j. Since the prefactors in (69) are indepen- 
dent of initial position, each x™J '* would simultaneously 
set the jth term in the power spectrum of every flux to 
its torus-averaged value. Evaluating the time-averaged 
fluxes for any of the individual initial conditions X^-j '* 
would not produce a torus-averaged flux. However, since 
the average of a sum of terms must equal the sum of 
the individual averages of those terms, the sum of the 
resulting | C^- : a 1 2 (each evaluated at a possibly different 



Or, rather, six times, since each X™q* wom d likely also be dif- 



Xo-J '*) would yield the torus-averaged value of all fluxes 



simultaneously. Recalling that each integrand in (68) is 



different anyway, there is no further waste in evaluating 
each one using a different initial condition X*o^- '*■ 

It is important to note that we have simply made an 
existence argument for Xo-e > Xq-l ano - evei T Xo-j 
What those values actually are would vary from problem 
to problem, and finding them for the Teukolsky prob- 
lem may not be practical. The integrands in (68) are 



not especially analytically transparent, so it may be that 
they can only be determined by evaluating those inte- 
grands for several initial conditions xoi which would de- 
feat the purpose of invoking the mean-value theorem in 
the first place. Still, we believe the potential added sav- 
ings from knowing the Xaj' * merits exploring whether 
the Teukolsky calculation harbors some structure or sym- 
metries that would allow those initial conditions to be 
determined with little or no added expense. We turn to 
those additional potential savings now. 



B. Relative cost of time-averaged vs. 
torus-averaged functions on low-order resonant tori 



Assume that we have in hand the x .j ' for each j and 
agree to evaluate the coefficients Z* m -. x using those spe- 
cial initial conditions. The added efficiency is twofold: 
each Cj-x should potentially be less expensive to com- 
pute than any given A^ n (by reducing a double integral 
to a single integral), and fewer such Cj- X s than AknS 
will have to be computed in order to achieve a given tar- 
get accuracy in the torus-averaged fluxes (by reducing 
a double sum to a single sum). In fact, the more effi- 
cient calculation might even increase the resulting flux 
accuracy. We justify those claims in turn below. 



1. Cost of a coefficient 

For ease of illustration, we will estimate the relative 
computational costs of a single Cj ]X and of any single 
Akn for which k,n satisfy the selection rule (48). To 



make the comparison more stark, we remap the integral 
(68) to the interval [0, 2ir] via a linear change of variable 



Xp = OpA 



(72) 



to obtain 



7* = — 



2tv 



dXpe^fL r ,x(x(Xp;Xo)) ■ (73) 



The relative cost of the single integral ( 73 1 and its 



double-index counterpart (61 ) will depend on the specific 



ferent from 



and 



X 0:L Z 



numerical integration algorithms used to evaluate them 
and are difficult to estimate. However, we can sketch a 
crude argument that the single integral should be more 
cost efficient by considering the Fast Fourier transform 
(FFT) as the algorithm. 
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Consider first the ID integral (73 1, which we can inter- 



pret as the jth Fourier coefficient of a periodic function 
on [0, 2ir}. For a periodic function, an FFT will return all 
the Fourier coefficients from C-m 1 through Cn ± by sam- 
pling the integrand at 2iVi + l equally spaced abscissae 23 . 
So to capture Cj , we would need 2 | j \ + 1 evaluations of 
the integrand. However, the highest index coefficients 
computed via an FFT are heavily afflicted by aliasing 
error, while the lowest index coefficients computed are 
relatively free of such error. To minimize aliasing effects, 
we imagine increasing the number of sample points (and 
thus of coefficients computed) by some integer safety fac- 
tor 24 S so that Cj will be one of lowest index coefficients 
returned by the FFT and thus fairly free of aliasing er- 
ror. The total number of integrand evaluations under 
this scheme for computing Cj would thus be S(2 \ j\ + 1). 



Now imagine evaluating the double integral (61 1 using 



a 2D FFT, which we (even more crudely) envision simply 
as nested ID FFTs. Assuming the same safety factor 
S throughout, we would need S (2 \n\ + 1) <S (2 \k\ + 1) = 
S 2 (2 |n| + l)(2 I&I + 1) function evaluations. Re-expressing 
j in terms of n and k via the selection rule and using the 
number of integrand evaluations as a metric of numerical 



expense, the ratio of the cost of A kr 
would be 



to the cost of Ci 



cost of A kn (2|rc| + l)(2|fc| + l) 
cost of Cj 2 \nz + kp\ + 1 



(74) 



Generally speaking, for small values of both p and z, the 
denominator in the cost ratio is smaller than the numer- 
ator since n and k more often than not have opposite 
signs for a given j. It is conceivable that a single Ak n 
could turn out less costly to evaluate than Cj, but the 
likelihood of that would become higher as both p and z 
became large, for which case a resonant torus would be 
barely distinguishable from a non-resonant torus in terms 
of all the aspects discussed in this paper. 

The argument above artificially increases the true cost 
of evaluating both integrals and is not intended even to be 
fully convincing, let alone a proof. Rather, it is a heuristic 
illustration of a rule of thumb in numerical integration 
that, with similarly behaved integrands, ID integrals are 
less costly to compute than 2D integrals. 



2. Number of coefficients 

In contrast to the relative cost of computing a coeffi- 
cient, we can say more definitively that the total number 
of single-index coefficients needed to achieve some spec- 
ified accuracy in the torus-averaged fluxes will be less 



23 To make the formulae that follow more intelligible, we are sep- 
arately counting the value at 2ir, even though it is the same as 
the value at 

24 Ref. |31l recommends a factor of at least 4 for most applications. 



than the number of double-index coefficients needed to 
obtain the same accuracy. 

Suppose achieving a certain flux accuracy for a given 
I, m pair requires computing all A kn with indices up to 
n max and /c max - Denote the total number of torus co- 
efficients computed by J\[%. The Z* coefficients satisfy 



\ Z LJ = Zt(-m){-w) [22] ! so one of n and k need only 
run over non-negative values to obtain all the coefficients 
with |n| < n max , \k\ < /c max . The total number of coeffi- 
cients actually computed is therefore (having k run only 
non-negative) 



TVjj = (2n max + 1) (fc max + 1) - 



(75) 



For comparison, we determine the number Af\ of Cj 
coefficients (evaluated at x™J''*) that we would have to 
calcu late so that^ in light of the arguments of subsection 
VA every above would automatically be incl uded 

in the sum of all I (7,-r. As we showed in subsection V A 



the maximum j index that needs to be included in the 
single-index series that will thusly catch every k, n pair 



(76) 



The symmetry of Z* implies that j need not run both 
positive and negative, and the number of Cj 's we would 
need to calculate to ensure at least the same level of flux 
convergence as that attained with the Ak n coefficients is 



A/a 



vK 



l 



(77) 



Comparing equations ( 75 ) and ( 77 ) , we see that we need 
a factor of 



A4 



■A/a 



(78) 



(2& max -f" 1) (^max ~t~ 1) ^max 

^max ~t~ pkma,K ~~ 1 

fewer coefficients. The reduction in the number of co- 
efficients therefore depends on the order of the periodic 
orbit as well as on n max and fc max . The lower the values 
of p and z, the greater the reduction factor. 

Figure [7] showed the average number of kn modes on 
a resonant torus per distinct frequency. It also gives a 
general sense of how A4 av ings varies with fc max and n max . 
The agreement between the two is not exact because, 
when computing all j coefficients up to the maximum 
jmaxj some additional frequencies will be included that do 
not correspond to any of the included kn frequencies with 
|fc| < fc lriax , \n\ < n max . Therefore, Figure[7]overestimates 



but only slightly and gives a better estimate for 



the larger values of /c max , n m 
equation ( 78 ) gives A/" sa vin g s 
and A4avings ~ 12.45 for n, 



. For example, for q r g = |, 
7.84 for n max = fc max = 50 
ax = & max = 80, both of 



which agree with the values in the histogram of Figure [7] 
within a few percent. On the basis of Figure [7J we can 
therefore conclude that focusing on temporal rather than 
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spatial Fourier coefficients and invoking the above mean- 
value arguments could reduce by a factor of an order 
of magnitude or so the total number of Z* coefficients 
required to obtain accurate torus-averaged fluxes. 



VI. CONCLUSION 

Computation of adiabatic inspirals with a grid of reso- 
nant orbits could be an order of magnitude more efficient 
than the same computation with a non-resonant grid. 
If our speculations are verified and double sums can be 
collapsed to single sums (and double integrals to single 
integrals), there may be substantial additional savings 
since fewer and simpler coefficients will be required. To 
date, no accurate adiabatic EMRIs have been computed. 
Such a dramatic boost in speed would bring EMRIs more 
within computational reach. 
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Appendix A: Torus coordinates and time coordinates 

1. Mino time vs. coordinate time Fourier 
coefficients 

Suppose Akn-t is a Fourier coefficient of some bipcriodic 
coordinate time function f(t), 



A kn -.t = lim — 

T->oo T 



Then Ak n - t is also the knth Fourier coefficient of a dif- 
ferent biperiodic Mino time function 5(A), 



Akn-t = lim — 

A— >oo A 



(A2) 



.4 



kn: A 



We prove equation ( A2 ) by constructing the function g 
from / explicitly. We wil need the fact (see Ref. [55] for 
details) that dt/dX depends on r and 9 and is thus biperi- 
odic when evaluated on a trajectory r(X),0(X). dt/dX 
also has a nonzero average value T defined by 



r = lim — 

A-i-oo A 



A/2 dt 
dX 



-A/2 



dX 



(A3) 



Consequently, the function t(X) takes the form 



t(X) = TX + At(X) 



(A4) 



where Ai(A) is biperiodic in A and has zero average value. 

We can now construct g(X). Since Ai(A) is biperiodic, 
it is also bounded. Thus, if we define T = t(A), then in 
the limit T — > 00, we get T — > TA. We can now change 



(Al) variables in the integral in ( Al I 



T £A 

lim - I" dte l{k ^ +n ^ )t f{t)= lim \ [ ' dXe l ^ +nuj ^ rx+At ^ ^(X)f (t(X)) 

OO 1 A / TA aA 

J ~t (A5) 
= lim i [ 2 d Xe^ +nn ^ x e l ^ +n ^ At( - x ^(X)f(t(X)) . 

A-S-OO A I A dX 
2 



In the second line above, we have absorbed 25 the T into 
the A and used the fact that the coordinate-time and 



25 All we need is for the denominator of the prefactor and the size 
of the integration interval to agree. Since there is no preferred 
size for that interval (the function / is not periodic) and we are 
taking the infinite limit, we are free to call the size of that interval 
TA or A. 



Mino-time frequencies will be related [55] by 



UJg 



r 6 



(A6a) 
(A6b) 



to convert to Mino frequencies in the argument of the 
exponential. 



Comparing ( A2 ) to ( A5 1 , we see that due to the de- 



pendence on k and n in the argument of the exponential 
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(the coordinate time frequencies ujg and ui r appear here 
simply as parameters), there is actually a different func- 
tion 



3 fc „(A) = e^+^ Ai W^(A)/!/(A)) 



(A7) 



for each Akn:\- In other words, the Ak n -t are Fourier 
coefficients of a single function / while each Ak n - : \ is the 
fcnth Fourier coefficient of a different function gkn(X). 
But that poses no problem — we only sought to show that 
every i-Fourier coefficient is also the A-Fourier coefficient 
of some function of A. The pragmatic importance of 
this fact has to do with the evaluation of coefficients of 
torus functions, which, though stated in slightly different 
language, is the crux of the original argument in Ref. |26j 
and which we discuss in Section IA 21 

Analogous reasoning to the above in the resonant case 
with 



L0 r 



P 

Z 



(A8) 



and p, z relatively prime shows that each coefficient Cj-t 
of f(t) can likewise be considered a Mino time coefficient 
Cj-\ of some different function 



dt 



ffi(A)=e«^ A *W-(A)/(t(A)) 



(A9) 



The difference is that now each of Ai(A), dt/dX(X) and 
/ (t(A)) is a singly periodic function of A with period 



Ap = zA r = pAg = 



2tt 

Ulp 



(A10) 



Temporal Fourier coefficients can be calculated in the res- 
onant case without first having to convert to Mino time. 
We will nonetheless express and evaluate coefficients Cj ;t 
as coefficients Cj-\ in order to parallel the non-resonant 
case. 



respect to coordinate time t, there will be coordinates 
7 = (7,., 79 ) such that 



7 r (t) = u r t + -f rQ 

Jg(t) = LOgt + 7o 



(Alia) 
(Allb) 



Note that in any set of angle coordinates in which the 
trajectory velocities are constant with respect to some 
time parameter, orbit trajectories will all be lines with 
the same slope 1 + q r g. Such coordinate systems are 
nevertheless distinct: identical ordered pairs in two such 
coordinate systems will not, in general, correspond to the 
same point on the torus. 

Though not unique, the Xr~Xs coordinate system on 
Yg is nevertheless uniquely useful. Since each of 7 r ,7e 
would be a combination of \ r and xe > each point on the 
projected r-p r curve of an orbit would be labelled by a 
pair of values (7 r ,7e) rather than by a single value Xr, 
and likewise for the projected 9-pg curve. This mixing of 
radial and polar motions in each torus coordinate makes 
most calculations harder than they need to be, and the 
impetus behind Xr~Xe coordinates is precisely the conve- 
nience that flows from torus coordinates that separately 
shadow radial and polar motion. 

Still, we sometimes are interested in values of quan- 
tities averaged over the 7 coordinates. Luckily, by the 
correspondence between temporal Fourier coefficients of 
biperiodic functions and spatial Fourier coefficients of 



torus functions, equations (A2) and (A7) further estab- 
lish that if Ak n .^ is the fcnth coefficient of the torus 
function /(-f) associated with f(t), then it is also the 
fcnth coefficient Ak n -,x °f the torus function gkn{x) cor- 
responding to 3fe n (A). This is important for evaluat- 
ing torus coefficients in practice: it is usually very dif- 
ficult to go from a function / (r(t),p r (t), 9(t),pg(t)) to a 
form /(t) explicitly while it is straightforward to go from 
g(r(X),p r (X),8(X),p e (X)) to gffl. 



A-based vs. t-based torus coordinates 



Figure [5j represents T|. as a compact 27r-by-27r square 
in the Xr j Xe angle coordinates defined in Section II B 



Kerr geodesies trace out lines on this torus-square at con- 
stant velocity. With respect to any other time parameter, 
geodesic curves continue to be lines on the torus-square, 
but their parametric representations are not in general 
linear in time, nor are their velocities on the torus con- 
stant. 

For any choice of time parameter, however, there is 
always some set of coordinates on the torus such that 
geodesic motion on that torus is linear in that time pa- 
rameter and has constant velocity 26 . For instance, with 



Appendix B: A synopsis of the Teukolsky formalism 

Here we summarize some relevant aspects of the 
Teukolsky formalism as applied to the EMM problem. 
More details of this application can be found in numer- 
ous references, including [Ti l H8ti20 l 122 1 l30 l 133] . Our goal 
in this aapendix is to justify the expressi ons for the co- 



efficients zfj^ n and zf^°°^ in equations (|6l|) and (|68|), 



respectively. 



Darboux's theorem gaurantees that there is a way to write Hamil- 



ton's equations with respect to any evolution parameter. If the 
system is integrable, there will then exist a transformation to 
angle variables on the torus that increase linearly with respect 
to that evolution parameter 15 . 
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1. The Weyl scalar, ip4 



In 1972 Teukolsky derived the master equation [341 135] . 
a separable partial differential equation (PDE) whose so- 
lution describes the propagation in the Kerr spacetime 
of small perturbations to fields of different spin-weights 
s: scalar, electromagnetic and gravitational. Each so- 
lution to the master equation is a separable function 
which can be written as a multipole expansion. There 
are two computational approaches to solving the mas- 
ter equation for each multipole mode: the time-domain 
approach, which solves the resulting PDE directly, and 
the frequency-domain approach, which further Fourier 
expands the solutions. For the purposes of extracting flux 
information from gravitational perturbations, frequency- 
domain codes are the accuracy standard and the ones to 
which our savings proposal applies. We thus restrict our 
attention to the frequency-domain approaches to solving 
the master equation. 

The combined multipole-Fourier expanded perturba- 
tions take the form 

B t/}(t,r,6,tp) = 

/OO 
dojR lmuj (r) s SZ(d)e-^ t+mv , (Bl) 

.771 

where u> denotes the coordinate-time frequency of the 
perturbations at the field point due to the source. Each 
s ip is a function of the field point (t, r, 0, tp) at which 
we wish to evaluate the perturbation. The s marker in 
equation (Bl I is a "spin- weight parameter" |34) which de- 



notes the perturbation type. For gravitational radiation, 
s = —2, and = 4 where 



p = — (r — ia cos ( 



(B2) 



The functions i?; mw (r) and _2<S^(#) (described in the 
next subsections) each depend on the parameter u> as 
a consequence of the separation of variables procedure. 
When the source is a geodesic, w turns out to be a dis- 
crete variable composed of harmonics of the radial, polar 
and azimuthal frequencies of that geodesic. That discrete 
dependence can be expressed differently for non-resonant 
orbits, 



Wmfcn = mUJ v + kuj r + TlUlg 



and for resonant orbits, 



]L)p 



(B3a) 



(B3b) 



Once u) becomes a discrete variable, we can replace the 
integral over all possible u in equation (Bl I with a sum 



over either to, k, n or m, j for non-resonant and resonant 
sources, respectively. 

Because everything in this paper deals with gravita- 
tional spin-weighting, we henceforth omit all the —2 sub- 
scripts. The net result is that equation (Bl| becomes 



V^4 (t, r, 0, <p) = p 4 £ R lmUmkn (r) SZ mhn (0) e -*-»-*+*^ 

Imnk 

(B4a) 

for a non-resonant source and 

V>4 (t, r, 6, p>)=p'J2 R ^ mi (r) SZ m3 (0) e~^ t+ ^ 

Imj 

(B4b) 

for a resonant source. 

2. The Spheroidal Harmonics 

The functions SZ {&) with a spin- weight of s = — 2 are 
the gravitational (tensor) spheriodal harmonics, a gener- 
alization of the likewise spin-weighted spherical harmon- 
ics. These functions satisfy 35 



(acu) cos 2 + Aauj cos 9 — 
sin 9d0 \ d9 



Am cos + 4 



sin 



e, 



sz (0) 



= 



(B5) 



Cim are the eigenvalues for which equation ( B5 ) has 



solutions. Solving for SZ{8) for given l,m,u re- 
quires simultaneously determining an eigenvalue Ci m and 
the associated spheroidal harmonic. These eigenvalue- 
eigenfunction pairs can be computed in several different 
ways (see Refs. PUmiH]). 

The spheroidal harmonics satisfy several orthogonality 
relations. The one we will need in this paper is that, for 
fixed to and Li, 



SZ(0)SZ m (0)sm< 



1 

2V 



(B6) 



where the overbar denotes complex conjugation. We have 
chosen a normalization of as in Ref. [M] . 



3. The radial Teukolsky functions 

Solving for the radial functions Rimu (?*) is more diffi- 
cult. Ri muJ (r) satisfy the inhomogeneous radial Teukol- 
sky equation [35] 

T imu (r) = A 2 |- (r) ) ~ V lmu (r) R lmu (r) . 

(B7) 

The potential Vi mu (r) depends in part on the eigenvalue 
Cim of SZ{6), so equation (B5I must be solved before 



the homogeneous or inhomogeneous version of equation 
(B7) can be. 
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The source term Timui is built by, among other things, 
evaluating S f^ (9 ) and two homogeneous solutions 27 



R 



in/up 



(r) to (B7) along the geodesic source. Two general 

in/up 



methods are described for constructing (t). One 

approach integrates the homogeneous Teukolsky equa- 
tion (or, equivalently, the better numerically behaved 
Sasaki-Nakamura equation [35]) outward from the hori- 
zon. The other expands -R]^f" P (r) in terms of hyperge- 
ometric functions (see |23j and references therein) and 

evaluates R^^ (r) directly at certain points, possibly 
extrapolating its values to nearby points with series ex- 
pansions. Both approaches are fairly computationally 
costly. For detailed explanations on these different ap- 
proaches and how various numerical problems are cir- 
cumvented, see I35rl4"0] . We elaborate a bit more on 
the structure of the source term below. 



4. The quantities ZfJ^ 

With the homogeneous radial solutions in hand, 
the inhomogeneous Teukolsky equations can be solved 
using the method of variation of parameters 28 |4"T] . 
The radial functions Rimu (r) can be written as 

[HE! 131 [33 IH S3] 

R lmu (r) = Z[? n Jr)RZJr) + ^(r)i^(r) , (B8) 
where ZP^ (r) and Zf° (r) are defined as 



■ (r) 



dr' 



, Rr m Jr')T lm Ur') 



J lmuj 



(r) 



dr 



(B9) 



The constant c is related to the Wronskian of -Rj^ (r) 
and R^ u (r). r + is the larger root of A and is the radial 
coordinate of the black hole horizon (the smaller root is 
denoted r_). 

Because we are only interested in the radiation going 
into the black hole and being carried away to infinity, 
we are only concerned about the asymptotic behavior 
of Rimuj (r) as r — > r + and r — > oo. In fact, the ho- 
mogeneous basis solutions have been chosen to have the 
simplifying feature that 



Zwnu (r - 


> r+) 


= z l 


Zhnui( r ~ 


>r + ) 


= 


Zl%iu)( r ~ 


-> oo) 


= 


Zhnuj(. r ~ 


-> oo) 


= z x 



Imu) 



(BIO) 



Other basis solutions to the homogeneous equation exist, e.g. 
the out/down basis -R°„ t 1 f d ° Wn (r). For a summary, see [3] and 
references therein. 
28 Most references use the method of Green functions, but variation 
of parameters works as well. 



Note that we have used the same notation for the func- 
tions zfj^(r) and for the constants zfj^ representing 
their asymptotic values at oo and r + , respectively. As 
mentioned in Section II V Al the literature seems stuck 



with the rather backward notational convention that 
Zwnui i s nonvanishing at oo while Zfi^ is nonvanishing 
at r+. 

The radial functions as r — > oo and r — > r + thus be- 
come 



R Wn,u = Rlrnu(r -> Oo) 



7, H r 3 p iulr * 



oo) 



Rj**,,., — Ri' 



>7CO Din / 

Zoo a 2 „— ikr 



(Blla) 



(Bllb) 



where and k = uj — ma/(2r + ) and r* is the Kerr tortoise 
coordinate defined by 



* / \ 2r+ r - r 4 

r (r ) = r H in 

r + — r_ 2 

r" + a 2 



2r_ r — r_ 
In — - — 



dr 



The coefficients ZfJ^ are found from 



Z 



H/oo 
Imuj 



dr 



,RTril p ir')T lmu {r') 



(B12) 



(B13) 



The source function Timu {f) in (B13) is an integral of 
the form 



TlmLj(r) 



dt dnB mu {t,r,e,ip)S^e iu>t e- im ^ . (B14) 



The source term is derived in [44j , but we have written it 
borrowing the notation B mLJ from [T5j. All that matters 
for our purposes is that B mu is built out of a series of 
operations on the null tetrad components of the energy- 
momentum tensor of the orbiting particle and thus con- 
tains delta functions and derivatives of delta functions 
centered on the source geodesic. Thus, the dil integral 
can be evaluated, resulting in every 9 and tp in (B14) be- 



ing replaced with the source trajectories 9 s (t), ip s (t) (the 
subscript s here denotes "source" as opposed to a spin- 
weight as earlier in this Appendix). 

Delta functions S(r — r s (t)) and derivatives thereof 
still remain in (B14), along with the integration over t. 
When we plug (B14) into equation (B13l, we can switch 
the order of integration for r and t and use those re- 
maining delta functions in r (we rename r' to r now, for 
simplicity) to replace every r with r s (t). The net result 
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is that zfj^ takes the form 



7 H/oo 
J lmuj 



dte i U t e -im 9t { t)l H/oo {th6s {t)) 



(B15) 



The functions depend on r and 9 both directly and 

via a combination of elementary functions, the spheroidal 
harmonics, the homogeneous radial Tcukolsky functions, 
and various derivatives thereof. Explicit expressions can 
be found in several sources (see, for instance, [HI |2"3"] ). 

We will now use this form for Z?J°° to define the co- 
efficients zfH^ and Z^2°°\ to which the efficiency argu- 



Imkn 



ments of Sections IV an- 



and [V] apply, 



respectively. 



5. The quantities ZfJ^ n 



We now show that, when the source is a non-resonant 



orbit, the w-dependence of zfj^ takes the form 

ryH/cO \ rrH I CO <r , \ 

^Imu - 2^ ^Imkn-X \ U ~ ^mkn) , 
kn 



(B16) 



where 



^mkr 



are the coordinate-time harmonic frequen- 



cies defined in equation (B3a|. Note that even though 



equation (B16) has the form of the Fourier transform of 



an almost-periodic function with respect to coordinate 
time t, we are free to interpret the coefficients of the 
delta functions either as Fourier coefficients of a coordi- 
nate time function or as Fourier coefficients of a (differ- 
ent) Mino time function. 

For the reasons stated in Appendix [3] we opt for the 
latter and begin by rewriting ( B15[ ) as an integral over 
Mino time. Treating all source coordinates as functions 
of A and using equation (A4), we get 



r H/oo 
oo 



) 

(B17) 

We now use the fact (see Ref. [35]) that, like dt/dX, 
dip/dX depends on r and 9 and is thus biperiodic when 
evaluated on a trajectory r(A), 9(X). dip/dX has a nonzero 
average value Q v defined by 



£l v = lim i / dX ~ 

A-s-oo A J_A/2 dX 



(B18) 



Consequently, the function ip s (X) takes the form 

ip s {X) = fl v X + Aip s (X) (B19) 

where Aip s (X) is biperiodic in A and has zero average 
value. Like its radial and polar counterparts, fl v is re- 
lated to the coordinate-time frequency uj v via 



1 



(B20) 



More generally, coordinate-time and Mino-time frequen- 
cies are related by 



uj = „ Q 

r 



(B21) 



In light of equations ( |B19| -( [B2l] ), equation ( |B17| ) be- 



comes 



7 H/c 



dX 



W e-" nA ^l^(r s (X),e s (X)) 



(B22) 



Note that the coordinate-time frequency ui still appears 
as a pai 

iwAt(A) 



as a parameter in both and in the argument of 



Like At(A) and A<p(A), the functions T-fJ^ are biperi- 
odic in A. We define the biperiodic function 

/*£°(A) ee e^'We-^'W^ (r.(A)A(A)) 

(B23) 

and Fourier expand it as 

fSTW = E zfJZn^-^^" ■ (B24) 

k.n 

Inserting ( |B24 ) into equation ( B22 1 yields 



„h/oo _ sr^ yH/oc / ,> m\ -i{ m n v +kn e +nn r )\ 

Irau / j lrnuj;kn-X I 

k,n ' J -°° 

k,n 

= E Z im / ^o W „;fcn;A 27r<5 ( W ~ ^kn) ■ 
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The multipole index m and the Fourier indices k,n all 
do double duty by helping to specify the value of the 

Consequently, each 



parameter ui in ZfJ^L 



Z wi^=Lo mkn ;kn-\ is full y specified by the four integers 



Imuj— uj m kn \kn\X ' 



l^m^k^n and we can define 



fjH I 'oo _ 2 r?H /oo 
lmkn:X lmuj—Lo m k n ;kn;X 



(B26) 



H/oo 



We will absorb the factor of 27r into the function f ( 
in the integral defining zfj^ n . x . 

Since each ZfJ^° n .^ is a temporal Fourier coefficient, 
then by the arguments of Section |lIID 1[ we determine it 
by instead computing the corresponding spatial Fourier 
coefficient 



z H/oc 
Irakri 

1 

W) 2 Jo 



2 " d Xr d Xe e ik ™ e-*" f?J~ Umkn (Xr , Xe ) 



(B27) 
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of the associated torus function in the x torus coordi- 
nates. Despite our earlier notation, the function ifj^ 
actually depends not only on r and 8 but also on their 
conjugate momenta (since it depends on the r and ve- 
locities via the energy-momentum tensor of the particle). 

fimuj inherits this dependence, and it is thus appropri- 
ate to write it as a function of the torus coordinates x 
that need not have any special symmetries on the torus. 

As we show in Appendix [Cj the averaged fluxes re- 
quired to evaluate the RHS of the adiabatic equations 



for £ depend on the ZfJ^ defined in (|B27|). Thus, equa- 



tion ( B27 1 verifies equation (61), on which the savings 
arguments of Section [IV B| are based. 



The quantities Z^°° x 



We now show that when the source is a resonant orbit, 
the w-dependence of Zf*{°° instead takes the form 



7 H/oo 
J lmuj 



7 H/oo 



(B28) 



The relevant quantities, then, are those in equation (B30 ) 



with the parameter w in set to 

Paralleling the non-resonant case, each Z^H°° , , is 
fully specified by the three integers l,m,j. By absorbing 
the factor of 2ir into the functions 1^1°° in the integrand 



of (B30), we can define the notationally more compact 



coefficients 



7 H/oo _ 
J lmj;X 



7 H/oo 
J lmu}—ujr n 



(B32) 



By the construction above, each such Z^J°^ is given by 



Z 



H/oo 



dXe^f^fX;^) ■ (B33) 



As we show in Appendix [C] the time-averaged fluxes 
for £ from resonant orbits depend on the Z^°° x defined 



in (B33). Thus, equation (B33) verifies equation (68), on 



which the more speculative savings arguments of Section 
IVl are based. 



where w m j are now the coordinate-time harmonic fre- 



quencies defined in equation (B3b). As before, we are 



free to interpret the coefficients of the delta functions 
either as Fourier coefficients of a coordinate time func- 
tion or as Fourier coefficients of a (different) Mino time 
function and take the computationally tractable latter 
option. 



Equations (Bf7)-(B23) carry over to the resonant case 



with the difference that each of At(X), Aip(X),r(X),9(X) 
and thus is now singly periodic with period Ap. 

Thus ful^ can be Fourier expanded in harmonics of a 
single fundamental frequency ftp, 



/zm(T( A ;xo) 



EgH/oo 
lmuj;j;\ 



(Xo)e 



-ijflpX 



(B29) 



with 



1 



UA e Jlmu> 



(A;Xo) 



(B30) 

The functions fuj^ (X',Xo) are induced from some torus 
function, so by the arguments of Section |IHD 2| both 

they and the coefficients Z^I^^. x (xq) depend on initial 
positions, which we can represent compactly as a depen- 
dence on initial position xo orL the phase space torus. 



Inserting (B29) into the resonant version of equation 
(|B22|) yields 



3 
3 

E 

3 



yH/oo _ gH/oa 
Imu) / j lmu>;j;\ 

H/oo 

7 H/oo 
J lmcj— 
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(B3I) 
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Appendix C: Fluxes from the Teukolsky formalism 

In this appendix, we review how the apparatus of Ap- 
pendix IB] yields fluxes of conserved quantities. Several 
authors [3 [3 HI 1201 H2 US] implement this Fourier- 
domain formalism in TB codes to calculate the radiative 
£ fluxes at radial infinity and the horizon to determine 
how the inspiral evolves. We show how expressions for 
time-averaged (as opposed to torus-averaged) fluxes dif- 
fer between non-resonant and resonant orbits. 



1. Overview of flux calculation 

We will not refer to the Q flux, but restrict our dis- 
cussion to E and L z . To determine the evolution of an 
inspiral in orbital parameter space, we use the the E and 
L z fluxes at radial infinity and the horizon as proxies for 
the local self-force. This subsection gives an overview of 
the calculation for finding these fluxes. 

From the Weyl scalar ip4, the gravitational waveform 
and the E, L z and Q radiation fluxes can be calculated. 
Specifically, we can calculate the polarizations h + and 
h x of the metric perturbations at infinity (i.e. the GWs) 



ip4 



2~W 



(h + + ih x ) 



(CI) 



After integrating equation (CI ) twice to get h + +ih x (the 



integration constants are set to zero), we can calculate an 
effective GW stress-energy tensor at infinity [25] as 



rpOO 

1 a/3 — 



1 / dh + dh + dh x dh x \ 

16n\dx a dxP 8x a Ox? / several' 1 ' 



wavelengths 
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The average over several wavelengths signifies the follow- 
ing [46] . The stress-energy tensor of GWs contributes to 
curvature in the analogus way that the stress-energy ten- 
sor of matter does. However, since GWs are fluctuations 
in the metric itself, we cannot define their stress-energy 
tensor at a field point because, with only information at 
one point, we cannot distinguish between the curvature of 
the background spacetime and the contributions to that 
curvature from fluctuations on the spacetime. To dis- 
tinguish between the background and the effect of the 
fluctuations, there needs to be either a length or fre- 
quency scale separation between the two. In our case, as 
r — > oo, the background curvature scale is much greater 
than the wavelength of the fluctuations. Therefore, we 
average over several gravitational wavelengths in order 
to smooth out the fluctuations and determine only the 
secular contribution of such GWs to the curvature. 

While the background curvature scale is much greater 
than the GW wavelength, GW detectors look not for spa- 
tial fluctuations to the metric but rather temporal ones. 
Therefore, in the spirit of applicability to actual exper- 
iments, it is more useful to distinguish between a back- 
ground frequency (i.e. a reciprocal of a curvature scale 
in the timelike direction) and the much larger frequency 
of the fluctuations. Therefore, rather then average over 
several wavelengths, we average over several periods to 
isolate the net effect of the GWs [4"B] . 

The energy and angular momentum fluxes carried to 
radial infinity by GWs are related to the components of 

T? B by, 



/ 1 



T, 



dE 

dtdA 
dL z 

dtdA 



(C3) 



The E and L z fluxes are calculated by integrating equa- 
tion (C3) over a 2-sphere of radius r on a constant t 



space 



be calculated by 

dE\°° 
dt/ t 

dt 



lim — 

T^oo T 



lim — 



dE 
~dt 

dLz 
dt 



dt 



dt 



(C5) 



Alternatively, we can calculate the torus-averaged 
fluxes by taking the average of the time-averages for all 
geodesies with a given set of orbital parameters over all 
possible initial conditions, 



dE\°° _ 1 

~dt/ i ~~ ~(2n) 2 Jo 

Lz\°° 1 
dt 



2w r2TT 

dlr„ I d-yg 
Jo 

2w p2ir 

, , d lro i d le 
(2ir) Jo Jo 



dE 



V 

dt/t 

dLz 
dt 



(C6) 



An analougus procedure can be performed to calculate 
the fluxes at the horizon [351 HZ] ■ 



2. Fluxes from ip4, 



Combining equations (B4| with (Bll) we find that 



,/,oo _ -4 yH 3 iu mkn r* qaui mkrl /n\ -iui mkn t+imip 



Iraki 



(C7a) 



for non-resonant orbits, and 



Imj 



-iLj m jt-\-imt£> 



ike hypersurface 



(07b) 



f ) ^ / ^ 

- / ^* 



(04) 



The time-averaged fluxes from a given geodesic can then 



Similar expressions can be found for ip^ , but for brevity, 
we will only proceed with the detailed computations for 
r — > oo. Also, we only work out the details for the E flux, 
but the L z flux follows exactly the same perscription. 

Following the prescription laid out in section |C 1| we 
find that 



'dEy" 
dt ) 4?r 



\lmuj I'm'uj' / 



several 
periods 

(C8) 



where denote the discrete variables that are two- indexed for resonant orbits and three- indexed for non- 
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resonant orbits. Performing the <p integration in equation 
(C8) yields 27r8 mm i. We thus set m = m! everywhere. 



For non-resonant orbits, equation ( C8 ) becomes 



dt 



(C9) 



I'm I 

Imknl'k'n' ^ mKn ^ mK ' n ' J I several 

periods 



We are interested in an infinite time-average of equation 
(C9). Therefore, we can drop the average over several 



periods because the time-averaging process will smooth 
out the fluxes so that the period averaging will have no 
further effect once we have time-averaged. The time- 
average of a function corresponds to the constant term 
in a Fourier expansion, so the argument of the exponen- 
tial in t will need to be zero. Therefore, performing the 
inifinite time-average yields the added conditions k = k' 
and n = n' . This equates the frequencies everywhere, 
including in the spheroidal harmonics. We can therefore 
now perform the 9 integration using (B6) and get 1 = 1'. 



The result is that the infinite time average in the non- 
resonant case is 



dE\°° _ 
dt/ t 



E 

Imkn 



1 



yH I 
ImkniX 



(CIO) 



In Section |IIID 1| we saw that infinite time averages 
over non-resonant orbits are the same as torus averages 
over non-resonant tori with respect to the correspond- 
ing torus coordinates. The time average over t on the 
LHS of (CIO) thus corresponds to a torus average over 



the 7 torus coordinates. We also saw that 



7 h 

lmkn-X 



I III h i 



where the Zf mkn 



are spatial Fourier coefficients 



(for fixed l,m) with respect to the x torus coordinates. 
Therefore, the torus averaged energy flux is 



dE\°° 



E 



Ini.kn ^mkn 



\z, 



H 

Imkr 
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This expression is true for all tori, as torus averages are 
insensitive to whether the orbits on that torus are reso- 
nant or non-resonant. 

Analogous arguments lead to the angular momentum 
flux at infinity. Similar arguments to those above then 
lead to the corresponding fluxes at the horizon. The up- 
shot is that all the torus-averaged fluxes are given by 



dE\ H/co 

dLS H/ °° 
dt 



E 

Ira kn 



H/oc 
Imkr. 
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where af mkn 



E ^a; 3 
Imkn 
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J Imkn 



7 H/<x 
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(C12) 



1 and the details of cxf mkn can be found 
in reference [35]. We note that there is no residual de- 
pendence on the initial conditions xo- 



We return to equation ( C9 ) and evaluate it for a reso- 
nant orbit, 



'f )" -\ <E E ^^<..«-*-'' M '--' [ *-MC-«K- ■ (ci3) 

\lmj I 1 j' j j u j SC vcral 

periods 



As was the case with the non-resonant infinite time aver- 
age, we can drop the averaging over several periods. Ad- 
ditionally, the resonant time-average picks out the con- 
stant term in the Fourier expansion, which results when 
j = j' . Therefore, 



dE\ x 
dt/ t 



Imj ""V 
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found similarly and are, 
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The rest of the E and L z time-averaged fluxes can be We remark that unlike the torus-averaged fluxes, the 
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time-averaged fluxes of resonant orbits clearly depend 

on the initial conditions of the orbit, since as we saw 
2 



B6 



7 H/oo 
J lmj\\ 



is not the same for all initial 



in Section 
conditions. 

Alternatively, we can write the time-averaged fluxes of 
equation ( C15 ) explicitly in terms of the torus coefficients 



Z hiZ- From section 



HID 2 



Z, 
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we know that 
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Therefore, we can rewrite equation (C15) as 



irp\ H/oo H/oo 

f } = E E £^<^^e^("^')^ + ( fc - fe ') x -> (C17) 

' * Imkn k'n': mkn 
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The explicit initial condition dependence is now ev- 
ident. Notice that if we average the time-averaged flux 
expressions of equation ( Cf 7 ) over all possible initial con- 



ditions, we reproduce the torus-averaged fluxes of equa- 
tion 0111) 
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